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

Конвективные течения.... Вып. 2

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

И.И. Вертгейм1, Б.И. Мывникова1, М.А. Кумачков2

1 Институт механики сплошных сред УрО РАН, 614013, Пермь, Академика Королёва, 1

2 Пермский государственный университет,

614990, Пермь, ул. Букирева, 15

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

ВВЕДЕНИЕ

Как известно, термокапиллярная или термоконцентрационная конвекция в тонком слое жидкости характеризуется не только существованием течения вдоль свободной поверхности, но и деформацией самой поверхности [1]. Это явление привлекает пристальное внимание исследователей, а также имеет многочисленные тех-

© Вертгейм И.И., Мызникова Б.И., Кумачков М. А., 2005

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

Как свидетельствуют результаты экспериментальных и теоретических исследований, возникновение термокапиллярного течения и деформация границы могут быть вызваны облучением поверхности раздела. В [5-7] представлены экспериментальные данные изучения подобных эффектов вследствие воздействия сосредоточенного источника тепла, созданного излучением лазера или некогерентного светового излучения. В зависимости от параметров источника излучения и от стадии развития процесса в экспериментах наблюдались различные типы деформации границы раздела - вблизи горячего пятна выпуклый или вогнутый мениск, а также более сложные конфигурации поверхности с переменной кривизной. Структура конвективного течения, как правило, была одновихревой с радиальным растеканием жидкости от нагретой области, но при некоторых значениях параметров на периферии течения возникали вторичные вихри.

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

[8] в длинноволновом приближении, обобщено в [9] на случай горизонтально-неоднородного нагрева, вызванного локализованным источником тепла, а в [7] - с дополнительным учетом вертикальной неоднородности теплопотока.

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

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

Рассмотрим тонкий горизонтальный слой вязкой несжимаемой жидкости глубиной Ь (рис. 1).

Нижняя граница слоя 2 = 0 считается теплоизолированной, а граница 2 = Ь + И (х, у, ґ) - свободной и деформируемой. Предполагается, что в жидкости установлено распределение температуры, соответствующее однородному подогреву снизу. Это поле искажено наложением слабонеоднородного теплового потока.

Рис. 1. Г еометрия задачи

Система управляющих уравнений в терминах горизонтальной и (х, у, 2, ґ) и вертикальной ^ (х, у, 2, ґ) компонент скорости и температуры Т ( х, у, 2, ґ) жидкости имеет вид

Р (Э ,й + (и V) и + -щи'} = — Ур + А и + и" (1.1)

Р (Э^ + (иV)w + ww'} = — р' + Аw + ^"- О (1.2)

V- и + w' = 0, Э Т + и Ут + wҐ = АТ + Т" (1.3)

Уравнения (1.1) - (1.3) записаны в безразмерных переменных с помощью следующих единиц измерения: длины - Ь, времени -Ь2/с, скорости - с/Ь , давления - щ/Ь, температуры - АЬ , где А - абсолютная величина равновесного градиента температуры. Символы ?і,у, С, о обозначают коэффициенты динамической вязкости, кинематической вязкости, температуропроводности и поверхностного натяжения жидкости соответственно. Штрих означает дифференцирование по вертикальной координате 2 , а операторы V и А - по горизонтальным координатам х и у .

Чтобы сформулировать граничные условия для системы (1.1) -(1.3), представим функцию, характеризующую мощность источника тепла, в виде Q (х, у, z) = q (х, у ) • qz (z). Тогда получим:

На поверхности z =1 + ^х, у) выполнение условия баланса нормальных и тангенциальных компонент плотности потока импульса П ^ = ~P ^ (VI ^ + Vki) приводит к соотношениям:

В построенной модели присутствуют следующие безразмерные

ную тепловую неоднородность, причем Z = 0 означает, что qz (z) = const.

2. ДЛИННОВОЛНОВАЯ АППРОКСИМАЦИЯ

Будем предполагать, что горизонтальный масштаб термокапиллярного течения, поля температуры и деформации поверхности намного больше, чем вертикальный. Это допущение позволяет применить метод асимптотических разложений в приближении длинных волн, аналогично приведенному в [8-10]. Согласно этому подходу, переменные, описывающие возмущения основного состояния, раскладываются в ряды по малому параметру e, характеризующе-

z = 0: u = w = 0; T' = — 1 + q(x,y)-qz (0)

z =1 + h (x,y): ht + u Vh = w

(1.4)

T' =Vh VT — 1+q (x y)-qz (1 + h); qz (1) °C

Vh

\

V1 + (Vh)2 у

(Сa — MaT)

v

(

\

(1.5)

1

параметры: число Марангони

Прандтля P = v /c, число Галилея

число

и капиллярное

число Ca = sLj(hc). Параметр Z в (1.4) характеризует вертикаль-

му отношение вертикального и горизонтального масштабов: г2 = (Э// Эх)/(Э//Эг) □ 1. В рассматриваемом случае для возмущений скорости, температуры, давления и деформации поверхности можно записать следующие разложения:

w = г(+ гм'1 + г2w2 +...), и =-1г\й0 + ги + г2и +...І (2.1)

Т = Т0 + е Т + е2 Т2 +..., h = е (Н + еМу +...) (2.2)

р = е 1Р-1 + Р0 + еР1 +... (2.3)

Проведем преобразование масштабов по формулам:

( х,у ) ^•\/е( х, у), / ®e2t, О ®еО, Са ®е2Са

В низших порядках по е, соответственно в минус первом и нулевом, из (2.1) - (2.3) получаем:

Р_1 = О(1 -z), Р0 = ОН - СаАН, Т0 = -z + Ф (2.4)

щ = 12/'(2)УФ, ^0 = -12/(z) АФ, /(z) = z2 - z3 (2.4)

Учет в разложениях (2.1) - (2.3) членов до второго порядка по е позволяет получить используемую в дальнейшем математическую модель задачи в длинноволновом приближении в виде системы нелинейных дифференциальных уравнений для функций Ф(х, у), ¥(х,у) и Н(х,у), характеризующих изменения амплитуд полей температуры, завихренности и деформации свободной поверхности вблизи критического значения числа Марангони Ыас = 48 :

Э,Ф + УФУх(е ¥) + — У4Ф-У2И-48у(|УФ|2Уф) + ' 1 15 35 ' >

- У (дУФ) +1 У2д +СчИ = 0

31?

е V2¥= —У(У2Ф)ХУФ-24УНхУФ

2 35Р 1 ’

ОУ2Н - СаV4Н = - 72У2Ф

(2.6)

(2.7)

Здесь в2 = (0,0,1).

Система (2.5) - (2.7) является обобщением модели, полученной в

[9], допускающим наличие вертикальной неоднородности локализованного теплового потока.

3. ОПИСАНИЕ НЕОДНОРОДНОСТИ ТЕПЛОВОГО ПОЛЯ

Предположим, что неоднородный источник тепла локализован в горизонтальной плоскости и расположен посередине слоя. В дальнейшем будут рассмотрены два варианта конфигурации теплового пятна. Поэтому функция q (х, у), описывающая горизонтальную неоднородность нагрева, будет зависеть от одной пространственной переменной, обозначаемой X и имеющей разный смысл в зависимости от изучаемой ситуации. Так что q = q (X), где X ° х в случае плоского источника и X ° г - для осесимметричного источника. Вычисления проведены для разных способов задания функции локальной тепловой неоднородности, показанных на рис. 2. Эти зависимости имеют вид гладких функций следующих двух типов:

Параметры функций q(X) характеризуют следующие свойства тепловой неоднородности в горизонтальном направлении:

- меру отклонения теплового потока от критического значения, соответствующего порогу возбуждения термокапиллярной конвекции, как внутри теплового пятна (- а2), так и вне его (р2);

где А =(а2 -р2^(Я/р))/(1 + ЩЯ/г)), В = р2 + А

(3.1)

(3.2)

- характерный горизонтальный размер тепловой неоднородности

я;

- ширину переходного слоя р по температуре между областями внутри теплового пятна и вне его.

Рис. 2. Варианты горизонтальной неоднородности теплопотока д(Х) при различных параметрах а, Д, р

Выбор вида функции д(X) в случае I связан с удобством сопоставления полученных результатов с ранее выполненным линейным анализом задачи [9], где предполагалась ступенчатая форма тепловой неоднородности. Ступенчатая конфигурация является предельным случаем зависимости I при р ® 0 . Вариант II, описываемый выражением (3.2), позволяет получить решение системы уравнений (2.5) - (2.7), которое в предельном случае недеформи-руемой поверхности сопоставляется с точным стационарным решением [10] и используется для тестирования численного метода.

4. РАСЧЕТЫ СТАЦИОНАРНЫХ РЕШЕНИЙ

Будем считать, что существует и является единственным стационарное решение системы (2.5) - (2.7), которое зависит только от координаты X и локализовано в пространстве так же, как и тепло-

вая неоднородность. Строгое доказательство этого факта в настоящее время отсутствует, хотя его достоверность косвенно подтверждается результатами многочисленных расчетов. Такое решение имеет вид:

Условие обращения в (4.1) завихренности в нуль следует из (2.5)

- (2.7) для рассматриваемых вариантов функции тепловой неоднородности.

Основное состояние (4.1) является решением системы уравнений

Использование переменной X подразумевает применение в (4.2)

- (4.3) преобразований на основе следующих тождеств: для плоского случая V/ ° (/',0) и У(g,0) ° g', для случая осесимметричного теплового пятна V(g,0) ° (rg)'/г . Здесь символы / и g обозначают любую скалярную функцию от X, штрих - производную по X, (Ц 0) в тождествах - вектор в ортогональной системе координат, где первая компонента направлена вдоль оси X.

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

численных расчетах рассматривается конечная область XI £ В ширины В □ К . Расчеты выполнялись с использованием метода мно-

(4.1)

ОУ2Н0 -СаУ4Н0 =-72У2Ф0

(4.3)

гократной пристрелки, реализованной в процедуре BVPMS библиотеки IMSL программного пакета Visual Fortran.

Рис. 3. Вид профилей Н0(£ для разных типов д(Х) при различных значениях параметра а(д= -0.015, 3 = 0.3): а - тип I, Р = 1, р = 0.1, 0.31 < а< 0.54; б - тип II, Р = 13, 0.7 < а< 5

Обсудим сначала поведение функций Ф0 (X), Н0(Х). На рис. 3 приведены зависимости Н0(Х) для плоского случая неоднородности, заданной выражениями (3.2), (4.1) при различных значениях параметра а и фиксированном значении параметров 3 и р . Расчет проводился в области шириной Б = 5 .

Видно, что существуют области изменения параметров а , 3, в которых происходит резкое изменение характера решения. Для них было выполнено более подробное исследование, результаты которого приведены на рис. 4 как зависимость Нт(а), где Нт °Н0(0)

- величина деформации поверхности в центре нагретой области.

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

Рассмотрим влияние параметра р для тепловой неоднородности типа II на вид решения. Результаты расчетов приведены на рис. 5. При относительно малых значениях р (тепловое пятно с резкой границей) решение соответствует подъему поверхности в центре пятна. Далее с возрастанием р наблюдается появление углубления на поверхности.

а б

Рис. 4. Изменение величины деформации поверхности в центре слоя для разных типов #(£) при различных значениях параметра а(д = -0.015, р = 0.3): а - тип I, Р = 1, р = 0.1, 0.01 £ а< 0.3; б - тип II, Р = 13, 0.7 £ а< 5

0.8

0.6

0.4

0.2

0

-0.2

-6

-4

-2

0

6

Рис. 5. Вид профилей Н0(Х) при различных значениях параметра р: тепловая неоднородность типа I, Р = 1, £=-0.015, а= 0.01, р = 0.01, 0.01 < р< 2.0 (кривые 1-9)

В рассматриваемой модели был введен параметр £, характеризующий изменение теплового потока по вертикали. В частности, случай д < 0 соответствует начальному этапу локального нагрева жидкости лазерным лучом, когда температура приповерхностного слоя жидкости выше, чем в глубине. При уменьшении |^| происходит выравнивание температуры по вертикали в процессе прогрева жидкости. Результаты, полученные при разных значениях данного параметра, представлены на рис. 6. Как видно, форма поверхности существенно меняется в зависимости от этого параметра задачи.

2

4

Часть результатов была получена с использованием квазиспек-трального метода расчета по полной двумерной системе уравнений

(2.5) - (2.7) [7]. Численные результаты показали, что найденные одномерные стационарные решения устойчивы и реализуются в определенном интервале изменения параметров. Качественно область устойчивости может быть определена из результатов линейной теории в приближении нулевого стационарного решения [9], для более точного определения надо рассматривать полную задачу устойчивости решений (4.2), (4.3). При выходе за пределы области устойчивости в двумерных расчетах могут развиваться как локализованные возмущения другой симметрии (в плоском случае - с нечетным профилем зависимостей Ф(X), Н(X), в осесимметричном

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

а б

Рис. 6. Вид профилей Н0(Х) при различных значениях параметра д(Р = 1): а - тепловая неоднородность типа I, а =0.01, ^ =0.01, -0.5 < д< 1.5; б - неоднородность типа II, а= 1.5, Ь= 1, -0.015< д<0.2

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

Заключение. Предложена теоретическая модель термокапиллярной конвекции, вызванной локально-неоднородным тепловым источником, и на ее основе выполнены расчеты стационарных режимов. Показана зависимость формы деформации поверхности,

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

Работа выполнена при частичной поддержке администрации Пермской области и РФФИ (гранты РФФИ-Урал № 04-01-96029 и РФФИ НШ-1981.2003).

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

1. Гершуни Г.З., Жуховицкий Е.М. Конвективная устойчивость несжимаемой жидкости. М.: Наука, 1972. 392 с.

2. Boeck Th., Karcher C. Low-Prandtl-number Marangoni convection driven by localized heating on the free surface: results of threedimensional direct simulations // Interfacial Fluid Dynamics and Transport Processes / R. Narayanan, D. Schwabe (Eds). Lecture Notes in Physics. 2003. 628. P. 157-175.

3. Полежаев В.И., Белло М.С., Верезуб Н.А. и др. Конвективные процессы в невесомости. М.: Наука, 1991. 240 с.

4. Безуглый БА., Иванова Н.А., Зуева А.Ю. Термокапиллярная деформация тонкого слоя жидкости, вызванная пучком лазера // Прикладная механика и техническая физика. 2001. Т. 42. № 3. С. 130-134.

5. Казенин Д.А., Карлов С.П., Шитиков Е.С. Нестационарная деформация свободной поверхности жидкости при локальном воздействии лазерного излучения // Труды Первой Российской национальной конференции по теплообмену. М.: МЭИ, 1994. Т. 6. Двухфазные течения. С. 96-99.

6. Мизев А.И. Экспериментальное исследование термокапиллярной конвекции, индуцированной локальной температурной неоднородностью вблизи поверхности жидкости. 2. Источник тепла, индуцированный излучением // Прикладная механика и техническая физика. 2004. Т. 45. № 5. С. 102-108

7. Karlov S.P., Kazenin D.A., Myznikova B.I., Wertgeim 1.1. Experimental and numerical study of the Marangoni convection due to localized laser heating // J. Nonequilibrium Thermodynamics. 2005. (In press).

8. Golovin A.A., Nepomnyashchy A.A. & Pismen LM. Pattern formation in large-scale Marangoni convection with deformable interface // Physica D. 1995. V. 81. P. 117-147.

9. Вертгейм И.И., Мызникова Б.И. Устойчивость и структурa термокапиллярного течения в горизонтальном слое с локализованным нагревом // Гидродинамика / Перм. гос. ун-т. Пермь, 2002. Вып. 13. С. 39-55.

10.Любимов Д.В., Черепанов А.А. Устойчивость конвективного течения, вызванного неоднородным нагревом // Конвективные течения / Перм. гос. пед. ун-т. Пермь, 1991. С. 17-26

NUMERICAL STUDY OF STEADY STATES OF THERMOCAPILLARY CONVECTION FROM THE CONCENTRATED SOURCE OF THE HEAT

I.I. Wertgeim, B.I. Myznikova, M.A. Kumachkov

Abstract. The nonlinear regimes of thermocapillary convection which are induced in a horizontal liquid layer by a concentrated heat source are investigated numerically. The layer surface is deformable. The numerical simulation is performed on a mathematical model which involves nonlinear partial differential equations for two-dimensional amplitude functions associated with distributions of temperature, vorticity, and surface deformation.

In the long-wave approximation the model describes the contribution of thermo-capillary convection to the heat transfer as well as the degree of the interface deformation. The proposed model generalizes the existing one by taking into account the heating inhomogeneity.