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

ДИНАМИЧЕСКИЙ МЕТОД ОПРЕДЕЛЕНИЯ ТЕПЛОФИЗИЧЕСКИХ СВОЙСТВ ЖИДКОСТИ

Л.Б. Директор, В.М. Зайченко, В.В. Качалов, И.Л. Майков (Объединенный институт высоких температур РАН)

Введение

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

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

В этих методах при обработке осциллограмм колебаний капли, как правило, используется линейное решение уравнения колебаний в приближении капиллярных волн (пренебрежимо малое влияние гравитации) - формулы Рэлея для определения капиллярной постоянной [4] и соотношения Кельвина для определения коэффициента вязкости [5]. Это накладывает определенные ограничения на условия проведения измерений: необходимо обеспечить малый размер капли и незначительное начальное возмущение (амплитуда колебаний не должна превышать 1 % радиуса капли), что, в свою очередь, требует использования прецизионной аппаратуры для регистрации изображений капли.

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

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

В работе предлагается экспериментальная методика определения коэффициентов поверхностного натяжения и вязкости на основе использования характеристик осцилляций большой капли жидкости при ее падении из капилляра на плоскую поверхность (рис. 1).

Рис. 1. Схема эксперимента: 1 - микрометрический винт; 2 - поршень; 3 - исследуемая жидкость; 4 - капилляр; 5 - подложка

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

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

Экспериментальная установка (рис. 2) состоит из измерительной ячейки с дозатором жидкости и сменными капиллярами, системы подсветки капли и скоростной цифровой видеокамеры. Фиксируется процесс отрыва капли от капилляра, ее падение и затухающие колебания формы капли на плоской поверхности. В зависимости от размера капилляра радиус сферы, эквивалентной объему капли, составляет 1,5-2,5 мм, характерная частота колебаний 50-70 Гц. Очевидно, что при таких частотах для построения адекватной осциллограммы колебаний необходимо иметь 15-20 точек за один период колебаний, что соответствует частоте съемки 1000 кадров/с. В экспериментах использовались цифровая камера Phantom 1540, которая при частоте 1000 кадров/с обеспечивает разрешение 1024x1024 пикселей, и цифровая камера VS-FAST, которая при той же частоте съемки обеспечивает разрешение 1024x500 пикселей.

Рис. 2. Схема экспериментальной установки: 1 - цифровая скоростная видеокамера; 2 - объектив; 3 - измерительная ячейка; 4 - источник света; 5 - компьютер

Обработка результатов экспериментов

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

Исходная информация представляет графические файлы в стандартном точечном формате ВйМаР (ВМР), каждый из которых содержит изображение формы капли в разные моменты времени (рис. 3 а).

Рис. 3. Изображение капли: а - растровое изображение; б - контур меридионального сечения

На первом этапе обработки серое изображение переводится в черно-белое монохромное, при чем выбирается среднее значение цвета из палитры цветов. В результате определяется граница белого и черного цветов и, соответственно, координаты точек контура изображения (рис. 3б). Выбор граничного цвета при переводе изображения из серого в монохромное вносит определенную ошибку в результат. Так, например, при выборе пятой части полной палитры относительная погрешность определения объема капли составляет 0,2 %.

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

высокоскоростной цифровой видеокамеры фиксировались процессы отрыва капли от капилляра, падения капли и фаза затухающих колебаний капли на подложке (рис. 4). Регистрация изображений производилась с частотой 1000 Гц (1 мс соответствует одному кадру) на ПЗС-матрице размером 1024^1024 пикселей.

£ * • Л. V т •

0с 2х10"2 с 2,8x10'2 с 3,0х10"2 с 3,7хЮ'2 с 6,9x10'2 с

Рис. 4. Кадры характерных фаз колебаний капли воды

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

В результате обработки изображений определяются основные геометрические характеристики капли: радиус, площадь поверхности и объем капли как функции времени. Средний объем вычисляется осреднением соответствующей осциллограммы. Максимальная относительная ошибка определения объема в эксперименте составляет 0,2 %. Соответственно, ошибка в определении радиуса эквивалентной сферы составит для данного разрешения матрицы (масштабный коэффициент 0,0159 мм/пк) менее

0,1 %. Для расчета среднего периода колебаний использовалась сплайн-интерполяция [9].

С точки зрения уменьшения объема вычислительных процедур логично было бы период колебания и коэффициент затухания определять по осциллограмме колебаний вершины капли. Однако при большой начальной амплитуде и значительных размерах капли (приближение капиллярных волн неприменимо) колебания

со

О 200 400 600 800 1000 1200

/, мс

Рис. 5. Осциллограмма колебаний площади поверхности капли

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

Расчет коэффициента поверхностного натяжения

Рассмотрим каплю идеальной несжимаемой жидкости, лежащую на горизонтальной плоскости. Предположим, что жидкость не смачивает поверхность (угол смачивания равен п). Уравнение возмущенной поверхности осесимметричной капли в сферической системе координат (0, г) с началом в центре тяжести капли и осью г (осью симметрии), перпендикулярной горизонтальной плоскости, запишется в виде:

Г (0, ()=я0 (1+((0)+^(0,0), (1)

где Л(0) - не зависящая от времени функция, описывающая равновесную форму капли в поле силы тяжести; £(0, 0 - осесимметричное возмущение поверхности капли; Я0 - радиус сферы

эквивалентного объема капли. Будем считать, что |£(0, О|<<|Л(0)|. Форма капли описывается функцией

Р(г,0,^ = г - гР (0,^ = 0. (2)

Изменение поверхности капли во времени записывается в виде:

dt dt

(U V)F = 0,

(3)

где и - вектор скорости.

Уравнение движения капли имеет вид:

ти ^ /л,

р— = -УР -pge, (4)

т

где р - давление; р - плотность жидкости; g - ускорение свободного падения; е - единичный вектор в направлении оси г.

Из уравнений (3) и (4) можно получить соотношение

Э2 ^ Эр 1 д^ Эр

dt2 dr

Э0 Э0

(5)

где pt = ^ — Pn (cos 0)аи (t) - зависящая от времени составля-

п V R J

ющая давления; Pn(cos0) - полиномы Лежандра; an(t) - коэффициенты разложения в ряд.

Подставив выражение (2) в уравнение (5), получим уравнение:

д

( - Y

PR° dt2 ^R° I R°

R

Pnan

( - Y

R°

d(R + ^)dP„ dPn

Э0

Э0

an, (6)

справедливое на границе капли при г = Я0 (1 + Я + ^).

Умножим уравнение (6) на Рт 8Ш0 и проинтегрируем по углу от 0 до п:

Э 2

-Р^2 д1 (Щ1, П ) =

ЭЯ Э

= х(((п) + nt- 1)I(m,Я,n))-)^m,——,n \an, где ^ = ^Pn (cos 0) (t); в„ - коэффициенты разложения;

n

I(m, f,n) = JPJ()Pn sin0d0.

n

Граничное условие для давления на поверхности капли запишется в виде:

Ро -PgZ + Z П ( + R )Pn an =

_ n (8)

= r-( + (L - 2)R + (L -2)£-2^LR - 2RLE).

где z = rF cos 0 = R ( + R + ^)) (cos 0);Pj(cos0) - полином Лежандра 1-го порядка; о - коэффициент поверхностного натяжения жидкости, L =-----1—— sin 0— - угловая часть оператора Лапласа

sin 0 Э0 Э0 и LPn = n(n + 1) Pn.

Динамика поверхности соответствует нестационарной части уравнения (8):

-Р^с I PnPP +1 n (1 + R )Pn а и =

( " " N (9)

= RI ( - 2)ZPnPn - 2^enPnLR - 2RL^PnP ]■

-**0 V n n n J

Умножив уравнение (9) на Pm sin0 и проинтегрировав по углу от 0 до п, получим:

^aп (I(m,1,n) + nI(m,R,n))-

n

a f (I(m,1,n)-2I(m,R,n)- ^ (10)

в ^ (n - i)(n + 2) V У ’ = 0,

n R I I 2I(m,LR,n) + BoI(m,_P1;n)

PgR

где Во = - число Бонда.

а

Следует отметить, что суммирование в формулах (7) и (9) начинается с п = 2. Это связано с выполнением условий постоян-Г 4 з

ства объема капли J <ЗУ = 3 и неподвижности ее центра масс

| МУ = 0 [10].

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

,2 О °» = Т^

^ п(п -1)2п +1) ^

п +^1 {т, Я, п )-

РЯ3

2п +1 ( дЯ д

--------I1 т,--------, п

2 I д0д0

х

(п - 1)(п + 2 )(1 -(2п +1)7 (т, Я, п ))-

-(2п +1)7 (т, ЬЯ, п) + Во

х

1 + •

г(2п +1)

4п2 + 4п -1 4п2 + 4п - 3 /

Л-1

х

I (т, Я, п)

(11)

х

Минимально возможная частота колебаний капли соответствует п = 2. При Я = 0 и Во = 0 формула (11) переходит в формулу Релея [4]:

ю2 = -80-. (12)

Р%

Пренебрегая вкладом формы капли в уравнении (11), получим соотношение

ю2 = -80- (1 + Во—1, (13)

Р*31 84

которое описывает гравитационную поправку к частоте собственных колебаний.

Результаты обработки экспериментальных измерений и расчета коэффициента поверхностного натяжения при различных числах Во показали, что только при числах Во < 0,01 ошибка в определении коэффициента поверхностного натяжения по соотношению Рэлея составляет менее 1 %. Условие Во < 0,01 соответствует ^о<005 [2с

— < 0 , 05 , где а =---капиллярная постоянная.

а ]|pg

Введение поправки (13) позволяет вычислять коэффициент поверхностного натяжения с приемлемой точностью до значений числа Во = 0,35.

Расчет коэффициента вязкости

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

у=а + °2 ехР

V °3 у

(14)

где а1, а2, а3 - параметры; t - время.

Коэффициент вязкости определяется по формуле [5]:

Гр^о2

Г 1

где 1 =—.

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

(15)

где у. - координаты экспериментальных точек; у - координаты расчетных точек.

Расчетные точки у. являются функциями параметров а1, а2, а3:

у = а + а2 ехр

V аз у

(16)

Разложим (16) в ряд Тейлора в окрестности точки (0, а0):

Уі = У і ((, 0°, о°, Оз0) + А« + ехр

п £

V «з у

Ао2

(17)

2 „2

ехр

V «3 у

Аа3.

5

а

і=1

Минимум функционала (15) определяется из условий:

дЬ = о,

д(Аа1) дЬ

д(Аа2) '

дЬ д(Аа3)

= о, = 0.

(18)

Подставив (15) в (18) и продифференцировав, систему уравнений (18) можно записать в виде:

У« - У -А°1 - ехР

V °3 у

Аа2 +

о ^

+ а2-2ехр

а3

V 3

(

V °3 у

Аа3

= 0,

У« - Уг -Аа1 - еХР

V а3 у

Аа2 +

+ а0 —2 ехр

а3

V 3

(

V а3 У

Аа3

ехр

V а3 У

= 0,

У, - Уг -Аа1 - еХР

X,

V а3 У

Аа2 +

+ а2-2ехр

2 2 а3

V 3

V а3 У

Аа3

о

а2 —2 ехр

V а3 У

= 0,

1=1

г=1

г=1

или

3 С N Эу Л N

Е Е ^ Аак = Е У - Уг )

к=1 V г=1 дак У і=1

^ 12 ІТяр

к=1 ^ і=1 дак

N

= Е У - Уг )ХР

V аз у;

^ак =

і=1 3 ( N

\ аз У

1II-^а? А,ХР

к=1 ^ г=1 дак а3

= Е У - Уг К -ТЄХР

V аз уу

^ г. ^

^ак =

V аз У

(19)

Эу, , Эу,

где — = 1; — = ехр Эа1 да2

V аз У

Эу о (,

да3

■ - а2 — ехр

а.

V аз У

3 ^3

Новые значения а{0 (где к = 1-3) рассчитываются через найденные значения Аак по формуле

ак = ак +Аак ■

(20)

Программная реализация

Разработанные алгоритмы реализованы в программном комплексе, блок-схема которого представлена на рис. 6. Программный комплекс разработан в виде приложения для операционных систем Windows 98, Windows ХР с использованием среды программирования ББЬРШ. По функциональности программный комплекс можно разделить на три блока: блок обработки изображений, вычислительный блок и блок формирования результатов вычислений.

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

г=1

В результате этого этапа обработки находится контур капли в виде функции диаметра горизонтального сечения от высоты капли и записывается в текстовый файл. Время работы этого блока определяет время работы всей программы. Вычислительный блок состоит из программы обработки изображения эталона, расчета масштабного коэффициента и программ вычисления теплофизических свойств жидкости. В программе реализованы линейные модели расчета коэффициентов поверхностного натяжения (с учетом гравитационной поправки) и вязкости [5]. Возможно применение более сложных нелинейных моделей, но в этом случае значительно увеличивается время счета. Блок вывода результатов формирует файл отчета обработки экспериментальных данных.

Считывание графических файлов

Преобразование изображений в черно-белые монохромные

Считывание изображения эталона

Расчет масштабных коэффициентов

I

Сканирование пакета изображений капли

к

Определение радиуса эквивалентной сферы

1 ' 1

Расчет коэффициента Расчет коэффициента

поверхностного натяжения вязкости

1 1 г

Формирование протокола результатов обработки

Рис. 6. Блок-схема программного комплекса

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

Результаты

Тестовые эксперименты проводились с дистиллированной водой в атмосфере воздуха при комнатной температуре. В экспериментах варьировались диаметр капилляра и высота падения капли. Продолжительность процесса регистрации до момента времени, когда амплитуда колебаний уменьшалась до 1-2 пикселей, составляла 1,5-2 с.

Суммарная погрешность определения объема капли, амплитуды и частоты колебаний определяется разрешением ПЗС-матрицы, погрешностью измерения эталона, оптическими искажениями объектива камеры, точностью юстировки оптической схемы. При использовании цифровой камеры Phantom 1540 максимальная методическая погрешность определения коэффициента поверхностного натяжения составляет 2-3 %, коэффициента динамической вязкости 8-10 %.

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

Коэффициент поверхностного натяжения дистиллированной воды

Фаза эксперимента Способ обработки Коэффициент поверхностного натяжения, Н/м Отклонение от справочного значения, %

Летящая капля Формула Рэлея 0,0709 2,87

Колеблющаяся на подложке капля Формула Рэлея 0,0792 8,5

Формула Рэлея с поправкой на гравитацию 0,0729 0,14

Лежащая капля Уравнение Юнга-Лапласа 0,0702 3,82

Используемое в экспериментах с капиллярными волнами выражение для коэффициента вязкости опирается на аналитическое решение линейного уравнения колебаний (решение Кельвина [5]). Неучет нелинейных эффектов приводит к большой погрешности, а поправка на влияние гравитационных сил только незначительно корректирует линейные члены. Имеющиеся работы (например, [13]) определяют поправки в случае малых колебаний (малое отклонение от шарообразной формы), и величины поправок составляют доли процента. Разумные результаты получаются только на хвосте осциллограммы, где колебания приближаются к капиллярным. Расчет по нелинейной модели [2] дает существенно лучшие результаты. В различных сериях экспериментов значения коэффициентов динамической вязкости для воды лежало в диапазоне 0,9-0,95 мПас, т.е. отличие от справочного значения не превышало 10 %.

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

где ю - круговая частота колебаний; 1/х - коэффициент затухания колебаний.

При регистрации изображения с помощью цифровой техники амплитуда колебаний разрешается в пределах 1-2 пикселей, соответственно, условие разрешимости амплитуды последнего из фиксируемых колебаний

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

Если использовать известные решения для капиллярных волн:

Заключение

юх > 1,

юх > 2%п,

условие применимости можно записать как

Безразмерный комплекс в левой части неравенства является комбинацией критериев Re и We (We0,5/Re) и отражает соотношение вязких сил и сил поверхностного натяжения. В частности, при диаметре эквивалентной сферы 4 мм этот комплекс составляет: для воды - 0,0025, для ртути - 0,0004, для глицерина - 3,5. Соответственно для воды разрешается порядка 40 периодов колебаний, для ртути - 250, для глицерина динамический метод не применим - колебания апериодические.

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

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

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

1. Rhim Won-Kyu. Thermophiscal properties measurement of molten silicon by high-temperature electrostatic levitator: density, volume expansion, specific heat capacity, emissivity, surface tension and viscosity / Won-Kyu Rhim, Kenichi Ohsaka // Journal of Crystal Growth. - 2000. - Vol. 208. - P. 313.

2. Glorieux B. Surface Tension of Liquid Aluminia from Contactless Techniques / B. Glorieux, F. Millot, J.C. Rifflet // Int. Journal of Thermophisics. - 2002. - Vol. 23. - No 5. - P. 1249.

3. Meier W. Surface tension and viscosity of surfactant from the resonance of an oscillating drop / W. Meier, G. Greune, A. Meyboom, K.P. Hofmann // Eur. Biophis. J. - 2000. - No 29. - P. 113.

4. Ландау Л.Д. Гидродинамика / Л.Д. Ландау, Е.М. Лифшиц. -М.: Наука, 1986. - 736 с.

5. Ламб Г. Гидродинамика: пер. с 6-го англ. изд. / Г. Ламб. - М.: ОГИЗ-Гостехиздат, 1947. - 928 с.

6. Майков И.Л. Численное решение задачи о затухающих нелинейных колебаниях капли вязкой жидкости / И.Л. Майков, Л.Б. Директор // ЖЭТФ. - 2008. - Т. 133. - № 6. - С. 1314.

7. Майков И.Л. Численная модель колебаний капли жидкости / И. Л. Майков, Л.Б. Директор // Вычислительные методы и программирование. - 2009. - Т. 10. - С. 148-157.

8. Директор Л.Б. Численное моделирование динамики вязкой жидкости / Л.Б. Директор, И.Л. Майков // Изв. РАН. Механика жидкости и газа. - 2009. - № 5. - С. 100-108.

9. Форсайт Дж. Машинные методы математических вычислений / Дж. Форсайт, М. Малькольм, К. Моулер. - М.: Мир, 1980. -280 с.

10. Григорьев А.И. О некоторых свойствах разложений по производным от полиномов Лежандра, проявляющихся при исследовании нелинейных осцилляций капли вязкой жидкости / А.И. Григорьев, С.О. Ширяева // ЖТФ. - 2005. - 75 (9). - С. 20.

11. Директор Л.Б. Усовершенствованный метод лежащей капли для определения поверхностного натяжения жидкостей / Л.Б. Директор, В.М. Зайченко, И.Л. Майков. // ТВТ. - 2010. -Т. 48. - № 2. - С. 193-197.

12. Директор Л.Б. Экспериментальное определение коэффициентов динамической вязкости и поверхностного натяжения жидкости методом затухающих колебаний капли с большой начальной амплитудой / Л.Б. Директор, И.Л. Майков, А.А. Середа // В сб.: Теплофизические свойства веществ (жидкие металлы, сплавы и наносистемы); тр. II Международного семинара. - Нальчик: КБГУ, 2006. - С. 67-70.

13. Bratz A. Surface oscillations of electromagnetically levitated viscous metal droplets / A. Bratz, I. Egry // J. Fluid Mech. - 1995. -Vol. 298. - Р 341-359.