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

Физика

УДК 532.5 001: 10.14529/ттрИ160107

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

В.Г. Речкалов

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

Ключевые слова: моделирование; метод конечных элементов; капля; жидкость; испарение; термокапиллярная конвекция.

Введение

В работе [1] отмечается, что одним из существенных недостатков метода конечных разностей в задаче о свободном испарении капли жидкости является необходимость сглаживающей аппроксимации температуры вдоль свободной границы капли поскольку выбор аппроксимации содержит элемент субъективности. Вторым существенным недостатком является наличие «плохих» точек на свободной поверхности, расположенных слишком близко к регулярным узлам сетки, что может даже привести в некоторых случаях к расходимости процесса поиска стационарного решения. В частности нам так и не удалось получить решение для капли ацетона, а расчетные значения для скоростей точек потока в капле воды не соответствовали экспериментальным результатам. Все это породило обоснованные сомнения в правильности полученных решений.

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

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

Математическая формулировка задачи

Форма капли определяется из решения уравнения Лапласа для свободной поверхности жидкости. Движение жидкости в капле описывается уравнением Навье-Стокса:

— + (У-У)У = —-егаар + кАУ + , (1)

Ъ у ' о

где t - время; к - коэффициент кинематической вязкости жидкости; fbp =—- температур-

Ъо

0о ЪТ

ный коэффициент плотности; Т - температура; Т0 - температура основания капли; АТ = Т — Т0 ; g - ускорение свободного падения; р - плотность жидкости.

Для жидкости принимаем условие несжимаемости: = 0 .

1 Работа выполнена при поддержке фонда РФФИ, грант 13-03-00918-а

2 Речкалов Виктор Григорьевич - кандидат педагогических наук, доцент, кафедра общей и теоретической физики, Южно-Уральский государственный университет, г. Челябинск, Российская Федерация.

E-mail: viktor-rechkalov@mail.ru

Уравнение теплового баланса записываем сразу для стационарного состояния:

1 — —ДТ-(V-У)Т = 0,

(2)

где X и С - соответственно теплопроводность и теплоемкость жидкости;

Для исключения давления из уравнения Навье-Стокса возьмем ротор от его левой правой частей:

—+ юdivV + ^ -У)ю-(ю-У^ = к Дю-рД х егаа Т. дt

(3)

где со = го№ .

Учитывая, что для несжимаемой жидкости divV = 0 получаем общее уравнение для ротора: дю

— + ^-У)ю-(ю-Уу = кДю-рДхgradТ .

(4)

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

1 (д2Т 1 дТ д2Т1

—Т + +—2 дх х дх дг2

РС

[ д 2ю 1 дю д 2юл

—Т + +—т

дх х дх дг

дТ дТ

--и--^ = 0 .

дх дг

-| и+ ^ I+ ^ + ёРр^ = 0,

дю

дх

дю

[ д 2ю 1 дю д 2юл

—Т + +—т

дх х дх дг

дю

дх

дю

д

[ д2у +1 ду+д У

дх

2

дх дг

2

кю ию

хх

кю ию

хх

2 ду

-| ^ + ^ |+ _ + = —;

дТ дю дí

/

дх

хю = 0.

со -

ди дw дг дх

= 1ду = 1 ду

х дг х дх

(5)

(6)

(7)

(8) (9)

(10)

где и и w - проекции скорости; х и г - цилиндрические координаты; у - функция тока. Уравнение для ротора записано для стационарного (4) и нестационарного (5) случаев.

т-. т /- ию

В работе [2], которую мы упоминали, был потерян член — в уравнении для ротора.

х

Граничные условия

На рис. 1 показана область интегрирования с указание границ.

Функция тока у нас интересует с точностью до произвольной постоянной и, поскольку вдоль границы замкнутой области она сохраняет постоянное значение, удобно принять вдоль всей границы у = 0.

Ротор на оси симметрии (рис. 1, граница 3) равен нулю. Значения ротора на верхней границе

(рис. 1, граница 2) на каждом шаге итерации вычисляется по формуле (5).

На свободной поверхности капли (рис. 1, граница 1) ротор вычисляется по формуле [2]:

1 дадТ V адТ V ю =---+ 2— =--+ 2—,

тдТ дэ я т дэ я

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

Рис. 1. Область капли с указанием границ (повернуто). 1 - свободная поверхность, 2 - верхняя граница, 3 - граница вдоль оси капли

(11)

к

к

х

х

Речкалов В.Г. Моделирование конвекции в капле жидкости

при ее испарении методом конечных элементов

Температура на верхней границе (рис. 1, граница 2) считается заданной и постоянной.

На оси капли (рис. 1, граница 3)

— =0. (12)

Эх

На поверхности капли (рис. 1, граница 1) считается заданным тепловой поток

11ЭТпх +11ЭТпх =—Я , (13)

Эх Эг

где q - удельная плотность теплового потока на поверхности капли; пх и пг - направляющие векторы цилиндрической системы координат.

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

Схема алгоритма

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

Шаг 1. Нахождение температурного поля при значениях скоростей и и w, найденных на предыдущих итерациях (Уравнение 5). Вычисление ротора на поверхности капли (Уравнение 11).

Шаг 2, алгоритм 1. Нахождение ротора при известных значениях и и w и известных граничных значениях (Уравнение 6).

Шаг 2, алгоритм 2. Нахождение производной ротора с последующим вычислением приращений ротора за малый промежуток времени (ДО и новых значений ротора (Уравнение 7).

Шаг 3. Нахождение функции тока по найденным значениям ротора (Уравнение 8).

Шаг 4. Нахождение поля скоростей и новых граничных условий для ротора на верхней границе области (Уравнения 10 и 9).

Шаги 1-4 повторяются до тех пор, пока разность значений скорости в контрольной точке после двух последовательных итераций не станет меньше заданной величины. В качестве контрольной точки принималась точка на оси в средней части капли.

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

Алгоритм 1, являясь частной реализацией релаксационного метода [2], требует введения специального параметра или параметров для управления сходимостью. Пробные расчеты показали, что значения ротора на верхней границе капли претерпевают очень большие изменения, приводящие к беспорядочной смене направлений поиска решения. Для придания последовательности решений регулярного характера мы использовали «мягкое» изменение ротора на верхней границе области: щ = щ—1 + Лю , где с - новое значение ротора; щ—1 - значение ротора, полученное на предыдущей итерации; с - значение ротора, полученное из уравнения (9); число 0 < Л < 1, которое подбирается экспериментально, мы использовали в качестве необходимого параметра.

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

В большинстве случаев мы использовали оба алгоритма: на начальном этапе решения - алгоритм 1 и для уточнения решения - алгоритм 2. Использование обоих алгоритмов в паре повышает безотказность в работе программы и дает существенный выигрыш во времени.

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

Результаты моделирования

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

Ацетон

Ацетон

Бензол

Бензол

Толуол Рис. 2. Изотермы

Толуол Рис. 3. Функция тока

ИПС

ИПС

Вода

Вода

Ацетон

Бензол

ИПС

Вода

Толуол Рис. 4. Ротор скорости

1. Процессы, протекающие в капле жидкости при ее естественном испарении, только приближенно можно считать подобными. На рис. 2-4 можно наблюдать, как постепенно искажаются физические поля в соответствии с изменением скорости течения в капле. Особенно сильное ис-

Речкалов В.Г.

Моделирование конвекции в капле жидкости при ее испарении методом конечных элементов

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

2. Средняя скорость вихря на оси капли слабо зависит от размеров капли. Размеры капли нельзя изменять произвольно, соблюдая при этом геометрическое подобие. Чтобы получить капли разных размеров нам пришлось выполнить расчеты при различных значениях ускорения свободного падения. В табл. 1 приводятся значения для скорости течения при двух значениях ускорения: при естественном ускорении и при ускорении g = 5 м/с. При значении ускорения g = 5 м/с размеры капель увеличиваются примерно в 1,4 раза.

Таблица 1

Средние скорости течения (мм/с) на оси капли при различных значениях ускорения

g, м/с2 ИПС Вода Бензол Толуол Ацетон

9,8 2,80 3,05 13,49 14,14 118,30

5 2,95 3,23 15,40 16,20 122,50

3. Мощность поверхностных сил и сил вязкого трения составляет всего 10—5 % —10—4% от мощности теплового потока. Полная же кинетическая энергия вихря еще на два порядка меньше и составляет, например, для капли ИПС К = 1,5 • 10—8 мДж.

Приближенные формулы

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

1. Кинетическая энергия конвективного вихря внутри капли жидкости

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

K » —pR w , где R - радиус капли в средней ее части. Расчет по этой формуле для ИПС дает 16

результат К = 1,36 • 10—8 мДж, что достаточно близко к значению приведенному ранее.

2. Мощность поверхностных сил

^ов » X ^, (14)

где т - напряжения на свободной поверхности капли, v - скорость, Дя - площадь элемента поверхности.

йо ОГ ДГ

т =--»а—, (15)

йГ Д1

где Д1 - длина образующей элемента поверхности; ДГ - приращение температуры вдоль образующей.

Площадь элемента поверхности находится по формуле Дя = 2 ж х Д1 .

ДГ

Nnов » X а—2пхДЬ = 2аХ XVДГ . (16)

Заменим приближенно х » хср » ^, где R - радиус капли в ее средней части, и V » w , где w -

средняя скорость течения на оси капли.

Nnов » pRawX ДГ = pRawДГ . (17)

Для капли ИПС по результатам расчета программы мощность поверхностных сил равна Nnов = 2,66 • 10—6 мВт. Приближенная формула дает значение равное Nnов = 2,35 • 10—6 мВт .

В выражении справа А Т имеет смысл разности температур между основанием капли и ее вершиной. Считая процессы в каплях жидкости приближенно подобными, на основе анализа раз-

мерностей можно записать: ДТ

д

рм>С

1-

1

N

Н рСм, даЯ Г

, где Н - высота капли. Следовательно

1

\

рС у НрСм

3. Мощность сил вязкого трения можно найти по формуле[4]:

2 / л

N

тр

2 г Д^ 2

Дх

„и" А Дм

+2 и?+2 У дм

+

Ди Дм

— + —

Дг Дх

2

ду~и—V, ЯН

(18)

(19)

где V - объем капли.

4. Средняя скорость течения на оси капли

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

даЯ

(

1 --

1

2

-V.

(20)

рС V" НрСм) ' ЯН Из полученной формулы вытекает приближенное выражение для средней скорости вихря на оси капли:

V

да Я2 Н ( 1 - V д 1

трС V Н рСм )

(21)

где А - безразмерный коэффициент.

Безразмерный комплекс физических и геометрических параметров в круглых скобках предоставляет собой обратное значение числа Пекле: -= Ре- . Число Пекле характеризует соот-

Н рСм>

ношение между конвективным и молекулярным процессами переноса тепла. Для процессов, в которых конвективный перенос тепла преобладает над теплопроводностью число Пекле принимает большие значения. Для рассмотренных нами задач число Пекле принимает значения, лежащие в диапазоне от 64 для воды и до 4300 для ацетона (табл. 2). Следовательно, выражение в скобках можно считать приближенно равным единице.

Ограничивая рассмотрение только близкими по форме висячими каплями, формулу (21) можно упростить.

м » Б.

да

\ирС

где Б - безразмерный коэффициент.

Полученную формулу можно записать в другой форме —» Б-

д ¡¡рС

(22)

, в которой тепловая на-

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

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

Табл. 2 показывает, что коэффициенты А и Б действительно очень мало отличаются. По значениям коэффициентов жидкости можно разбить на три группы подобия: вода, ИПС; бензол, толуол; ацетон. Если посмотреть на рис. 2-3, то можно согласиться с тем что, процессы в каждой из групп максимально похожи. Наибольшие отличия мы наблюдаем в капле ацетона, значения коэффициентов которого сильно отличаются от значений в двух других группах.

Речкалов В.Г. Моделирование конвекции в капле жидкости

при ее испарении методом конечных элементов

Таблица 2

Значения коэффициентов А, В и числа Ре при заданных значениях q_

Жидкость д, мВт/мм2 А В Ре

Вода 0,10 1,781 1,730 75

ИПС 0,18 1,645 1,599 149

Бензол 0,25 2,407 2,353 496

Толуол 0,25 2,497 2,448 581

Ацетон 4,00 3,902 3,812 4207

Экспериментальная проверка

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

Рис. 2. Конвективный вихрь в капле ИПС. а) Фотография следов частиц, движущихся в капле жидкости. б) Линии тока, численный эксперимент

Рис. 2 демонстрирует принципиальное сходство экспериментальных и расчетных линий тока в капле изопропилового спирта. При этом надо иметь ввиду, что криволинейная поверхность капли создает оптическое искажение особенно сильное вблизи контура капли. Средняя скорость частиц, измеренная по результатам видеосъемки, составляет в среднем 2,5 мм/с. Расчетное значение средней скорости на оси капли равна 3 мм/с, что вполне согласуется с экспериментом.

Эксперименты с каплями воды дали неожиданный результат: никакого видимого упорядоченного движения обнаружить не удалось. Расчетное значение скорости движения в капле воды очень близко к скорости движения в капле ИПС. Мы ожидали увидеть примерно то же самое, что и на рис. 2, однако на фотографиях, которые мы не приводим, видимых следов упорядоченного движения нет. Приходится предположить, что на поверхности капли воды кроме термокапиллярных сил присутствуют иные силы, препятствующие развитию конвекции. Эти силы могут быть связаны с изменением свойств поверхности, вызванным поглощением ПАВ из окружающего каплю воздуха. Однако наличие ПАВ на поверхности капли должно было бы привести к изменению ее формы, что экспериментально не подтверждается.

Основные результаты

Получены уравнения термокапиллярной и гравитационной конвекции в капле жидкости при ее естественном испарении для конечных элементов.

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

Создана программа полу-автоматического построения однородной конечно-элементной сетки.

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

Выполнена экспериментальная проверка теоретической модели.

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

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

Обнаружено не объясненное пока отсутствие конвекционной активности в капле воды.

Литература

1. Речкалов, В.Г. Особенности моделирования процессов свободного испарения капли жидкости методом конечных разностей / В. Г. Речкалов // Наука ЮУрГУ: материалы 66-й научной конференции, секции естественных наук. - Издательский центр ЮУрГУ, 2014. - С. 267-274.

2. Бараш, Л. Ю. Испарение и динамика лежащей на подложке капли: дис. ... канд. физ.-мат. наук / Л. Ю. Бараш. - М., 2009. - 74 с.

3. Берковский Б. М. Вычислительный эксперимент в конвекции / Б.М. Берковский, В. К. Полевиков. - Минск.: Университетское, 1988. - 167 с.

4. Берд, Р. Явления переноса / Р. Берд, В. Стьюарт, Е. Лайтфут. - М.: Химия, 1974. - 688 с.

Поступила в редакцию 19 октября 2015 г.

Bulletin of the South Ural State University Series "Mathematics. Mechanics. Physics" _2016, vol. 8, no. 1, pp. 49-56

DOI: 10.14529/mmph160107

MODELING OF CONVECTION IN THE LIQUID DROP BY ITS EVAPORATION USING FINITE ELEMENT METHOD

V.G. Rechkalov1

The results of the modeling of thermocapillary and gravitational convection in a drop of liquid by finite element method are given. The approximate semi-empirical formula for a set of parameters of convective vortex are described. The experimental verification of the modeling results is made.

Keywords: modeling; finite element method; a drop of liquid; evaporation; thermocapillary convection.

References

1. Rechkalov V. G. Osobennosti modelirovaniya protsessov svobodnogo ispareniya kapli zhidkosti metodom konechnykh raznostey [Features of modeling of processes of free evaporation of a liquid drop by finite difference method]. Nauka YuUrGU: materialy 66-y nauchnoy konferentsii, sektsii estestven-nykh nauk [SUSU Science: Proceedings of the 66th Scientific Conference, Section of Natural Sciences]. Izdatel'skiy tsentr YuUrGU Publ., 2014, pp. 267-274. (in Russ.).

2. Barash L.Yu. Isparenie i dinamika lezhashchey na podlozhke kapli: dis. ... kand. f.-mat. nauk [Evaporation and dynamics underlying drops on the substrate. Cand. phys. and math. sci. diss.]. Moscow, 2009, 74 p. (in Russ.).

3. Berkovskiy B.M., Polevikov V.K. Vychislitel'nyy eksperiment v konvektsii [Computer experiment in convection]. Minsk, Universitetskoe Publ., 1988, 167 p. (in Russ.).

4. Berd R., St'yuart V., Laytfut E. Yavleniya perenosa [Transport Phenomena]. Moscow, Khimiya Publ., 1974, 688 p. (in Russ.). [Bird R.B., Stewart W.E., Lightfoot E.N. Transport Phenomena. New York, John Wiley and Sons, Inc., 1960, 780 p. DOI: 10.1002/aic.690070245]

Received October 19, 2015

1 Rechkalov Viktor Grigorevich is Cand. Sci. (Pedagogical), Associate Professor, General and Theoretical Physics Department, South Ural State University, Chelyabinsk, Russia. E-mail: viktor-rechkalov@mail.ru