ЖЭТФ, 2021, том 159, вып. 2, стр. 359-370

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

БЕНАРА - МАРАНГОНИ В ИСПАРЯЮЩИХСЯ ЛЕТУЧИХ

КАПЛЯХ НА НАГРЕТОЙ ПОДЛОЖКЕ

А. А. Гаврилина^a,b, Л. Ю. Барашb*

^a Национальный исследовательский университет «Высшая школа экономики»

101000, Москва, Россия

^b Институт теоретической физики им. Л. Д. Ландау Российской академии наук

142432, Черноголовка, Московская обл., Россия

Поступила в редакцию 24 августа 2020 г.,

после переработки 14 сентября 2020 г.

Принята к публикации 28 сентября 2020 г.

Изучаются нестационарные внутренние течения в лежащей на нагретой подложке капле капиллярного

размера, испаряющейся в режиме пиннинга контактной линии. Проведено трехмерное моделирование

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

потоков Марангони необходимо учитывать диффузию паров в воздухе, теплопередачу во всех трех фазах

и тепловое излучение. Уравнения были решены численно методом конечных элементов с использовани-

ем программного обеспечения ANSYS Fluent. В результате проведенных вычислений получено нестаци-

онарное поведение неустойчивостей Бенара - Марангони (БМ). На первом этапе появляется цветочная

структура ячеек БМ вблизи контактной линии. Для меньших контактных углов ячейки увеличиваются

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

ми экспериментальными и теоретическими исследованиями и помогают проанализировать и разрешить

связанные с этим вопросы.

DOI: 10.31857/S0044451021020152

капле из-за зависящего от температуры поверхност-

ного натяжения. Осесимметричные потоки Маран-

гони можно наблюдать при комнатной температуре

1. ВВЕДЕНИЕ

и атмосферном давлении [7-10].

Вопросы испарения капель привлекли значи-

тельное внимание не только с теоретической точ-

Термокапиллярный поток жидкости с нарушен-

ной осевой симметрией может самопроизвольно воз-

ки зрения, но и в связи с развитием новых прило-

жений, таких как образование структур на поверх-

никать внутри сильно нагретой капли с интенсив-

ным испарением и достаточно высокой скоростью

ностях, подготовка сверхчистых поверхностей, на-

несение микроматриц молекул ДНК и РНК, крис-

жидкости. В нагретой капле неустойчивость Ма-

рангони, возникающая из-за неоднородного распре-

таллография белков, сгорание капель топлива в

двигателе, диагностика заболеваний, разработка со-

деления температуры, может приводить к движе-

ниям жидкости, таким как продольные рулоны и

временных методов печати для струйных принте-

гидротермальные волны (ГТВ), что было впервые

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

продемонстрировано в [11] и подтверждено в мно-

устройств и др. [1-6].

гочисленных последующих исследованиях [12-19].

Неоднородный тепломассоперенос при испаре-

ГТВ — это распространяющиеся в жидкости перио-

нии капли имеет решающее влияние на кинетику ис-

дические во времени волны, вызванные неустойчи-

парения. В частности, он приводит к изменению рас-

востями в стационарной конвекции Марангони. ГТВ

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

были первоначально предсказаны для жидкого слоя

и может вызывать термокапиллярную конвекцию в

на твердой пластине Смитом и Дэвисом [20,21], а за-

^* E-mail: barash@itp.ac.ru

тем наблюдались в [22,23]. Согласно Смиту и Дэви-

359

А. А. Гаврилина, Л. Ю. Бараш

ЖЭТФ, том 159, вып. 2, 2021

су [20,21], возбуждение распространяется перпенди-

Сефиан и др. [16] продемонстрировали, что ГТВ

кулярно градиенту повехностной температуры при

в каплях FC-72 представляют собой волны, которые

малых числах Прандтля и в направлении против те-

распространяются по всему объему капли. Кроме

чения при больших числах Прандтля [20, 21]. Для

того, в работе [16] было изучено влияние ГТВ на

плоского слоя жидкости ГТВ и конвекция Бена-

теплопередачу в твердой подложке и на простран-

ра - Марангони (БМ) обычно происходят при каче-

ственное распределение скорости испарения. Осо-

ственно разных условиях. В то время как эффект

бый вид ГТВ был определен Карапецасом и др.

БМ становится более выраженным с увеличением

в [17], где наблюдались трехмерные спиралевидные

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

периодические во времени бегущие волны, организо-

пленках, появление ГТВ связано с горизонтальным

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

градиентом температуры (например, они могут воз-

углу. Ши и др. [24] определили критические числа

никать при боковом нагреве). В случае испарения

Марангони для возникновения неустойчивости Ма-

лежащих капель конкуренция между ГТВ и конвек-

рангони в лежащей на нагретой подложке капле с

цией БM более сложная.

низкой летучестью с помощью серии расчетов трех-

мерного компьютерного моделирования.

Первоначально ГТВ в каплях были идентифици-

Новые особенности рассматриваемой проблемы

рованы Сефианом и др. в 2008 г. [11] при исследова-

были выявлены Семеновым и др. [25], они обнару-

нии испарения лежащих капель этанола, метанола,

жили, что ячеистые структуры, найденные в кап-

воды и флюоринерта FC-72. Позднее были опубли-

ле с большим числом Прандтля, представляют со-

кованы более подробные наблюдения ГТВ для ка-

бой ячеистую структуру БM, а не ГТВ. Авторы на-

пель этанола на четырех разных подложках с раз-

блюдали неустойчивость Марангони в капле этано-

личной теплопроводностью [12]. ГТВ, движущиеся

ла капиллярного размера на нагретой подложке с

в азимутальном направлении, наблюдались для ле-

высокой теплопроводностью. Они также выполнили

тучих капель этанола и метанола. Для FC-72 кон-

численное моделирование при помощи односторон-

вективные ячейки возникали около вершины капли

ней модели, описанной в [26], где учтены, в частно-

и двигались к контактной линии. Было обнаружено,

сти, нестационарное испарение в диффузионной мо-

что количество волн возрастает с увеличением теп-

дели для неизотермической лежащей капли в режи-

лопроводности подложки и ее температуры и убыва-

ме пиннинга контактной линии и стефановское тече-

ет с уменьшением высоты капли. В то же время гид-

ние в газе. Структуры, полученные при помощи ре-

ротермальная неустойчивость в испаряющихся кап-

шения уравнения теплопроводности в капле и прове-

лях воды не была обнаружена.

дения численного моделирования, были в хорошем

В работе [13] в каплях метанола, этанола и FC-72

согласии с данными эксперимента для контактного

на подложке из политетрафторэтилена выделены

угла θ ≈ 29.3^◦. На первом этапе появлялся торо-

три стадии процесса испарения с термоконвектив-

идальный рулон, примыкающий к контактной ли-

ными неустойчивостями: нагрев капли (которая по-

нии; позже он дестабилизировался и распадался на

чти достигает температуры подложки), испарение

несколько трехмерных нестационарных ячеек БМ,

с термоконвективными неустойчивостями (тепловой

занимающих всю каплю. На третьей стадии ячейки

поток на этом этапе максимален) и испарение без

БМ сносились к контактной линии, где силы Ма-

термоконвективных структур (капля на этом этапе

рангони сильнее, и наблюдалась цветочная структу-

больше похожа на пленку). Основная фаза характе-

ра примерно с 15 ячейками БМ вдоль окружности,

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

что соответствует пространственному периоду, при-

ячеек связано с разностью температур между нагре-

мерно в два раза превышающему локальную тол-

той подложкой и температурой окружающей среды.

щину жидкости. Авторы обнаружили, что количе-

Временная эволюция ГТВ была также описана по-

ство ячеек БМ увеличивается с уменьшением высо-

средством одновременных наблюдений потоков теп-

ты капли.

ла и жидкости в испаряющейся капле этанола на

Взаимодействие конвективных ячеек БM и ГТВ

нагретой подложке [14]. Было продемонстрировано,

наблюдалось в недавней работе [27] в лежащей лету-

что ГТВ, распространяющиеся в азимутальном на-

чей капле силиконового масла на медной подложке

правлении в капле этанола, имеют место как в усло-

в режиме пиннинга контактной линии в широком

виях нормальной силы тяжести, так и в условиях

диапазоне температур подложки. Когда температу-

микрогравитации, так что влияние гравитации на

ра подложки была выше комнатной, наблюдалась

поведение ГТВ невелико [15].

квазистационарная конвективная ячеистая структу-

360

ЖЭТФ, том 159, вып. 2, 2021

Моделирование нестационарных неустойчивостей Бенара - Марангони. . .

ра БМ, переходящая в осциллирующее состояние

дель, которая учитывает теплопередачу во всех

по мере испарения. Когда подложка была холоднее

трех фазах, гидродинамику, диффузию паров в

окружающего воздуха, осциллирующие ячейки БМ

воздухе, стефановское течение в газе и тепловое

наблюдались, только когда температура подложки

излучение. С помощью программного комплекса

превышала 10.6^◦C. Размер ячеек становился боль-

конечно-элементного анализа ANSYS Fluent вы-

ше, а их количество уменьшалось, когда темпера-

полнено трехмерное компьютерное моделирование

тура подложки была выше температуры окружаю-

и получены нестационарные распределения тем-

щего газа. Для температуры подложки, превышаю-

пературы и поля скоростей в капле. В результате

щей 14.3^◦C, наблюдается сосуществование распро-

расчетов мы получили как цветочную структуру

страняющихся гидротермальных волн и колебатель-

ячеек БМ около контактной линии для достаточно

ной конвекции БM. В статье приведены критические

больших контактных углов, так и структуру, в ко-

значения краевого угла, при которых происходит пе-

торой более крупные ячейки БМ занимают границу

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

раздела жидкость-газ для меньших контактных

В частности, авторы [27] рассмотрели случай,

углов.

когда температура подложки на 7.81^◦C выше тем-

Статья организована следующим образом. В

пературы окружающей среды. В этом случае на пер-

разд. 2 описываются математическая модель и ме-

вом этапе появляется цветочная структура ячеек

тод численного решения. Значимость различных

БМ вблизи контактной линии. С уменьшением кон-

физических эффектов для рассматриваемых задач

тактного угла ячейки увеличивались в размерах и

также анализируется в разд. 2. В разд. 3 представ-

занимали центральную область поверхности капли

лены результаты и обсуждение. Раздел 4 содержит

(рис. 4 в [27]).

заключение.

Качественно похожие цветочные структуры яче-

ек БМ вблизи контактной линии на промежуточной

2. МАТЕМАТИЧЕСКАЯ МОДЕЛЬ

стадии испарения капли наблюдались в обеих пуб-

ликациях [25, 27]. В то же время пространственные

Рассмотрим трехмерную нестационарную зада-

структуры ячеек БМ при достаточно малых углах

чу для осесимметричной капли, лежащей на нагре-

смачивания, которые наблюдались в [27] и численно

той твердой подложке. Будем использовать прибли-

получены в [25], существенно различаются. Ожида-

жение поверхности капли сферической чашей. Это

ется сходство результатов, полученных в [25,27], по-

приближение допустимо при условии Bo ≪ 1, т. е.

скольку параметры капель, рассматриваемые в двух

когда влияние гравитации на форму поверхности

статьях, близки друг к другу. Также схожи реле-

мало. Здесь Bo = ρgh₀R/(2σ sin θ) — число Бонда,

вантные физические свойства 0.65 сСт силиконово-

h₀ — высота капли. Мы используем обозначения из

го масла и этанола, и в обеих работах использу-

таблицы.

ются нагретые подложки с высокой теплопроводно-

стью и близкими значениями температуры. Поэто-

2.1. Диффузия паров в воздухе

му удивительно, что результаты на поздней стадии

испарения капли сильно различаются. Целью насто-

Динамика концентрации пара u в окружающей

ящей работы является анализ поведения нестацио-

атмосфере описывается уравнением диффузии

нарных ячеек БМ в лежащих летучих каплях на на-

∂u

гретых подложках, включая стадию сравнительно

= DΔu.

(1)

∂t

малых углов смачивания. Также представляет инте-

рес определение основных физических механизмов,

Используются следующие граничные условия:

которые ответственны за исследуемое поведение и

u = u_s на поверхности капли, u = 0 вдали от капли,

должны учитываться при моделировании.

∂u/∂r = 0 и ∂u/∂z = 0 на осях соответственно r = 0

Мы теоретически исследовали структуры

и z = 0.

неустойчивости Марангони в испаряющихся ле-

Семенов и др. выводят в [26] и используют в [25]

жащих каплях этанола и силиконового масла

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

для параметров, соответствующих каждому из

дели нестационарного испарения лежащих капель,

экспериментов

[25, 27]. На основе трехмерного

которые согласуются с численными результатами

численного моделирования описаны потоки жид-

из [28]. Однако эти нестационарные эффекты ма-

кости и геометрические характеристики структур

лы на временных масштабах, превышающих время,

БМ в капле. Разработана математическая мо-

необходимое броуновской частице для прохождения

361

А. А. Гаврилина, Л. Ю. Бараш

ЖЭТФ, том 159, вып. 2, 2021

Таблица. Обозначения и значения параметров, используемые в вычислениях

Единица

0.65 сСт

Обозначение

Величина

Этанол

измерения

силиконовое масло

Радиус контактной линии

мм

2.95

2.09

Контактный угол

градусы

29.2^◦, 20^◦, 10^◦

28.39^◦, 20.49^◦, 17.18^◦

T_s

Температура подложки

307.05

302.69

T_a

Температура окружающего

297.55

294.88

воздуха

Плотность

кг/м³

772.24

760

c_p

Удельная теплоемкость

Дж/(кг·К)

2602.3

2000

Теплопроводность

Вт/(м·К)

0.14

0.1

κ = k/(ρc_p)

Температуропроводность

м²/с

6.97 · 10-8

6.58 · 10-8

k_air

Теплопроводность воздуха

Вт/(м·К)

0.026

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

кг/(м·с)

1.095 · 10-3

4.94 · 10-4

Поверхностное натяжение

мН/м

2.062 · 10-2

1.54 · 10-2

σ_T

Производная поверхностного

Н/(м·К)

-8.979 · 10-5

-8 · 10-5

натяжения по температуре

Удельная теплота

Дж/кг

918600

223000

парообразования

Коэффициент диффузии

м²/с

11.81 · 10-6

5.6 · 10-6

паров в воздухе

u_s

Плотность насыщенных паров

кг/м³

0.219664

0.469145

Коэффициент излучения

0.92

0.91

⎛

характерной длины R, т. е. при t ≳ R²/D ≈ 0.7 c, для

√

Du_s

J (r) = |D∇u| =

⎝sinθ+

2(x(r)+ cos θ)3/2

рассматриваемых задач. Поэтому кинетику испаре-

ния можно рассматривать как стационарный про-

⎞

цесс и использовать квазистационарное приближе-

∫

∞

ch (θτ)

ние

τ th((π - θ)τ)P-1/2+iτ (x(r)) dτ^⎠ ,

(3)

ch (πτ)

Δu = 0

(2)

где

√

r² cosθ/R² +

1 - r2 sin2 θ/R

с теми же граничными условиями вместе с прибли-

x(r) =

(4)

1 - r²/R²

жением сферической чаши.

Последняя задача математически эквивалентна

и P-1/2+iτ(x) — полином Лежандра. Выражение (3)

задаче об электростатическом потенциале заряжен-

удобно аппроксимировать как [30, 31]

ного проводника, форма которого определяется дву-

(

)-λ(θ)

мя пересекающимися сферами. Такая задача была

r²

J (r) = J₀(θ)

(5)

решена аналитически в [29] в тороидальных коорди-

R²

натах. На основе этого аналитическое решение урав-

где λ(θ) = 1/2-θ/π. Коэффициент J₀(θ) в (5) может

нения (2) с указанными выше граничными условия-

быть найден следующим образом:

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

тока испарения J(r) с поверхности испаряющейся

J₀(θ)

= J₀(π/2)(0.27θ² + 1.3),

(6)

капли, было представлено в работе [30] в виде

1 - Λ(θ)

362

ЖЭТФ, том 159, вып. 2, 2021

Моделирование нестационарных неустойчивостей Бенара - Марангони. . .

(

)

Λ(θ) = 0.22(θ - π/4)² + 0.36,

(7)

u_g

J (r) =

J (r),

(11)

u_s

u_g - u_s

Du_s

J₀(π/2) =

(8)

где u_g = M_gp_g/(¯T), pg — атмосферное давление,

где u_s — плотность насыщенного пара на свободной

M_g = M_air(1 - X_v) + M_vX_v — молярная масса га-

поверхности капли. В работе [31] численно найде-

за, находящегося рядом со свободной поверхностью

но, что уравнения (5)-(8) представляют собой весь-

капли, M_air и M_v — молярные массы воздуха и па-

ма точную аппроксимацию решения уравнения (2) с

ра, X_v = u_s

RT/(M_vp_g) — мольная доля пара. Здесь

указанными выше граничными условиями в прибли-

температура поверхности T является функцией r.

жении сферической чаши. Это было подтвержде-

Отсюда следует, что cтефановское течение в газе

но для широкого диапазона углов смачивания. При

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

этом точность аппроксимации не зависит от таких

и 18 % для капель этанола и 0.65 сСт силиконового

параметров, как D и R. По этой причине мы исполь-

масла соответственно.

зуем уравнения (5)-(8) при моделировании.

Результаты, приведенные ниже в разд. 3, не учи-

Давление насыщенного пара определяется урав-

тывают стефановское течение. Однако наши предва-

нением Антуана (с p_s [Па] и T [K]):

рительные численные расчеты для капли 0.65 сСт

силиконового масла с параметрами, указанными в

log₁₀ p_s (T ) = A -

(9)

таблице, и краевым углом 20.49^◦ показывают, что

C+T

влияние стефановского течения на распределение

где A, B, C — эмпирические константы. Мы исполь-

температуры на поверхности капли весьма мало.

зуем A = 10.247, B = 1599.039, C = -46.391 для

этанола [32] и A = 9.418, B = 1509, C = -31.55 для

2.3. Гидродинамика

0.65 сСт силиконового масла [18].

Концентрация пара u_s на поверхности капли рас-

Основными уравнениями гидродинамики внутри

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

капли являются уравнения Навье - Стокса и уравне-

ется по формуле

ние непрерывности для несжимаемой жидкости:

Mps (T)

∂v

u_s (T) =

(10)

+ (v · ∇)v +

grad p = νΔv,

(12)

∂t

где

R — универсальная газовая постоянная, M —

div v = 0,

(13)

молярная масса пара.

Здесь ν = η/ρ — кинематическая вязкость, v — ско-

Из выражения (5) следует, что потеря массы

рость жидкости, p — давление.

в результате испарения происходит неравномерно

Граничные условия для уравнений (12),

(13)

вдоль свободной поверхности жидкого слоя, значи-

включают условие прилипания v = 0 на границе

тельно увеличиваясь вблизи закрепленной контакт-

раздела подложка-жидкость и условие на границе

ной линии. Неоднородный поток массы во время

раздела жидкость-газ, включающее силы Марано-

испарения и соответствующая теплопередача изме-

гони, связанные с температурной зависимостью по-

няют распределение температуры вдоль поверхно-

верхностного натяжения. Вообще говоря, последнее

сти капли и вызывают силы Марангони из-за зави-

граничное условие можно выразить как [34]

сящего от температуры поверхностного натяжения.

(

))

Эти силы приводят к термокапиллярной конвекции

p-p_v -σ

n_i =

внутри капли.

R₁

R₂

)

⁽∂v

∂v_k

∂σ

=η

n_k -

(14)

2.2. Стефановское течение

∂x_k

∂x_i

Поскольку концентрация насыщенного пара не

где единичный вектор n направлен в сторону пара

является величиной, намного меньшей плотности

вдоль нормали к поверхности. Взяв тангенциальную

воздуха (см. таблицу), естественно учесть в рассмат-

составляющую этого уравнения, находим

риваемой задаче вклад стефановского течения в га-

)

dσ

⁽∂v_τ

∂v_n

dϕ

зовой фазе [33].

=η

-v_τ

(15)

dτ

∂n

∂τ

dτ

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

скорость испарения увеличивается следующим об-

Здесь dσ/dτ = σ_T dT/dτ, τ_i — компоненты единич-

разом [26]:

ного вектора τ , тангенциального к поверхности, ϕ —

363

А. А. Гаврилина, Л. Ю. Бараш

ЖЭТФ, том 159, вып. 2, 2021

угол между нормалью к поверхности капли и вер-

в окружающем газе, k_air — теплопроводность воз-

тикальной осью. Отметим, что пакет моделирова-

духа, J (r) — локальная скорость испарения, опре-

ния ANSYS Fluent позволяет напрямую задавать

деляемая выражением (5),

напряжение Марангони на свободной поверхности

˙q_r = σ_Bε(T⁴ - T4a)

(20)

капли, поэтому уравнение (15) уже реализовано в

программном обеспечении, и пользователю не нуж-

— поток излучения между каплей и окружающей

но его задавать.

средой, σ_B — постоянная Стефана - Больцмана, ε —

коэффициент излучения, T и T_a — температура по-

верхности капли и температура окружающего газа

2.4. Теплопередача в капле

соответственно.

Распределение температуры вычисляется с по-

Соотношение (20) можно использовать при сле-

мощью уравнения теплопроводности

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

мер капли был намного больше характерных длин

∂T

+ v · ∇T = κΔT.

(16)

волн теплового излучения, которые определяются

∂t

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

тельная способность практически не зависит от час-

Теплопроводность внутри подложки определяет-

тоты [35]. И наоборот, для объектов, размеры кото-

ся уравнением

∂T

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

= κΔT.

(17)

∂t

можна суперпланковская теплопередача излучени-

ем в дальней зоне [36,37]. Во-вторых, важно, чтобы

Граничные условия для уравнения теплопровод-

изменение температуры внутри капли было намного

ности включают изотермическое граничное усло-

меньше, чем разница температур капли и окружаю-

вие на нижней границе подложки T

= T_s, адиа-

щего газа.

батическое граничное условие на границе раздела

Получим теперь соотношение для (∂T/∂n)_air.

подложка-газ ∂T_s/∂n = 0 и непрерывность тепло-

Квазистационарное распределение температуры

вого потока на границе раздела подложка-жидкость

в газовой фазе определяется уравнением

k_s∂T_s/∂z = k∂T_L/∂z. Здесь T_L и T_s соответствуют

температуре жидкости и температуре подложки, k

ΔT = 0

(21)

и k_s — теплопроводности жидкости и подложки со-

с граничными условиями: T = T_s на поверхности

ответственно.

капли (в предположении, что разница температур

Согласно экспериментальным условиям [25, 27],

внутри капли намного меньше, чем разница темпе-

в расчетах в настоящей работе используется под-

ратур капли и окружающего газа), T = T_a вдали

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

от капли. Эта задача математически эквивалент-

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

на уравнению (2) с соответствующими граничны-

T = T_s на границе раздела подложка-жидкость.

ми условиями, поэтому можно использовать те же

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

формулы для аналитического решения. Из (5)-(8)

ности капли имеет вид

следует, что величина J(r)/(Du_s) зависит только от

∂T

Q₀(r)

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

(18)

∂n

ности следует, что

)

где

J (r)

⁽∂T

(22)

⁽∂T⁾

Du_s

T_s - T_a

∂n_air

Q₀(r) = -k_air

J(r) +

q_r

(19)

∂n_air

Следовательно, основное граничное условие

представляет собой полный тепловой поток на по-

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

(19), где

верхности капли, который включает скрытую теп-

(∂T/∂n)_air вычисляется при помощи формул,

лоту парообразования, эффекты теплопроводности

аналогичных (5)-(8):

в газе и тепловое излучение. Приведенные ниже

(

)-λ(θ)

⁽∂T⁾

r²

оценки показывают, что все эти эффекты важны

= J₁(θ)

(23)

∂n_air

R²

для капель, испаряющихся с нагретых подложек.

Здесь n — вектор нормали к свободной поверхно-

(

)

J₁(θ)

сти капли, k_air

∂T/∂n _air — тепловой поток, связан-

= J₁(π/2)(0.27θ² + 1.3),

(24)

ный с неоднородным распределением температуры

1 - Λ(θ)

364

ЖЭТФ, том 159, вып. 2, 2021

Моделирование нестационарных неустойчивостей Бенара - Марангони. . .

T -T_a

2.7. Численное моделирование

J₁(π/2) =

(25)

Мы решаем уравнение теплопроводности и урав-

нения Навье - Стокса методом конечных элементов

2.5. Теплопередача в газовой фазе

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

Оценим относительную силу тепловых потоков

ния гидродинамики ANSYS Fluent. Моделирование

LJ(r) и k_air (∂T/∂n)_air, которые являются слагае-

включает несколько шагов.

мыми в (19). Из (22) следует, что тепловые потоки

• Построение геометрии и сетки.

равны, когда

Создаем трехмерную геометрию осесиммет-

LDu_s

T_s - T_a =

(26)

ричной капли. Радиус кривизны поверхности и

k_air

высота капли равны соответственно R/ sin θ и

Следовательно, вклад эффекта теплопроводнос-

R(1/ sin θ - ctg θ). Следующим шагом является

ти в воздухе соответствует тепловому потоку, свя-

разделение области моделирования на малые

занному с испарением, для

вычислительные ячейки. Используем количество

ячеек от 200 до 500 тысяч в зависимости от угла

T_s - T_a ≈

смачивания.

⎧

• Определение модели и настройка параметров.

⎨

60.91 К для этанола,

Выбираем тип решателя со связью по давлению,

≈

0.94 К для гексанола,

(27)

⎩

нестационарный режим, вязкое (ламинарное) тече-

22.53 К для силиконового масла.

ние и включаем расчет энергии в модели. Указыва-

ем физические свойства жидкости: плотность, вяз-

В данной работе рассматривается случай T_s - T_a ≈

кость и температурную производную поверхностно-

≈ 10 К для капель этанола и силиконового масла,

го натяжения σ_T , позволяющую задавать напряже-

поэтому рассмотрение эффекта теплообмена в газо-

ния Марангони и т. д.

вой фазе целесообразно для получения количествен-

• Установка граничных условий и метода реше-

ных результатов.

ния.

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

2.6. Тепловое излучение

ствии с разд. 2. Скорость потери тепла (19) ука-

зывается с помощью пользовательской функции

Оценим относительную силу теплового потока

(UDF), написанной на языке программирования Си.

LJ(r) испарительного охлаждения и теплового по-

Используется алгоритм решения SIMPLE (Semi-

тока теплового излучения (20) между каплей и

Implicit Method for Pressure Linked Equations — по-

окружающей средой.

лунеявный метод для уравнений со связью по дав-

Из уравнений (5)-(8) следует, что тепловые по-

лению) метода Патанкара [38], доступный в пакете

токи примерно равны при

Fluent. Основной особенностью этого метода явля-

)1/4

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

⁽LDu_s

T ∼

+T4a

(28)

к скорости и давлению для обеспечения сохранения

σ_BεR

массы и вычисления поля давления. Алгоритм на-

Следовательно, вклад эффекта теплового излу-

писан таким образом, что уравнение неразрывности

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

выполняется автоматически.

ного охлаждения для

• Последующая обработка.

Мы отображаем результаты моделирования: по-

T_s - T_a ≈

ле скорости, абсолютные значения скорости и рас-

⎧

пределение температуры, и затем анализируем вих-

⎨

70.6 К для этанола,

ревую структуру и волновые картины.

≈

1.61 К для гексанола,

(29)

⎩

30.69 К для силиконового масла.

3. РЕЗУЛЬТАТЫ И ОБСУЖДЕНИЕ

В данной работе рассматривается случай T_s - T_a ≈

На рис. 1 показаны экспериментальные резуль-

≈ 10 K для капель этанола и силиконового масла,

таты из рис. 4 работы [27] и наши численные ре-

поэтому учет влияния теплового излучения целе-

зультаты для высыхающей капли 0.65 сСт силико-

сообразен для получения количественных результа-

нового масла. Здесь T_s = 302.69 К, T_a = 294.88 К,

тов.

R = 2.09 мм, а угол смачивания капли составляет

365

ЖЭТФ, том 159, вып. 2, 2021

Моделирование нестационарных неустойчивостей Бенара - Марангони. . .

Рис. 2. Временные зависимости неустойчивостей БМ в

капле силиконового масла для θ = 20.49^◦ . Графики со-

ответствуют следующим значениям времени испарения:

8 (а), 18 (б), 28 (в), 38 (г), 48 (д), 58 (е), 68 (ж) с. Вре-

менная эволюция профиля капли не учитывается, т. е. кон-

тактный угол остается постоянным во время моделирова-

ния

28.39^◦, 20.49^◦ и 17.18^◦. Диффузия паров в воздухе,

Из рис. 2 видно, что осцилляции структуры яче-

гидродинамика в капле, эффекты теплопроводности

ек являются очень медленными по сравнению с ха-

во всех трех фазах и тепловое излучение учитыва-

рактерными частотами ГТВ в [11, 12]. Таким обра-

лись в численных расчетах в соответствии с обсуж-

зом, полученные результаты подтверждают в дан-

дением в разд. 2.

ном случае наличие нерегулярной осциллирующей

конвекции БМ, а не ГТВ. Ячейки конвекции Маран-

Полученное поведение ячеек БМ качественно

гони являются достаточно глубокими и практически

согласуется с экспериментальными данными: при

достигают дна капли, как показано на рис. 3.

большом угле смачивания наблюдается цветочная

структура, а при малых контактных углах ячейки

Рисунок 4 был получен для испаряющейся капли

БМ занимают всю поверхность капли. При модели-

этанола. Здесь T_s = 307.05 К, T_a = 297.55 К, R =

ровании не учитывалась временная динамика фор-

= 2.95 мм, θ = 29.2^◦. Временная динамика профиля

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

капли не учитывалась. Численные результаты ка-

ственные результаты, такие как число ячеек БМ,

чественно согласуются с экспериментальными дан-

могут отличаться.

ными. Мы получили 20 ячеек БМ, что согласуется

Отметим, что число Прандтля довольно вели-

с численными результатами работы [25], получен-

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

ными другим методом. Мы подтверждаем резуль-

неустойчивости ГТВ [20, 21] имеет тенденцию пере-

тат [25] о том, что полученная неустойчивость долж-

мещать ячейки БМ в направлении от вершины кап-

на быть описана как нестационарная конвекция БМ,

ли к контактной линии.

а не как ГТВ.

367

А. А. Гаврилина, Л. Ю. Бараш

ЖЭТФ, том 159, вып. 2, 2021

Рис. 3. Вид сбоку для моделирования капли силиконового масла с θ = 20.49^◦

Рис. 4. Структуры в распределении температуры на поверхности и график векторного поля скорости для капли этанола

с θ = 29.2^◦. а) Экспериментальные данные, перепечатанные из работы [25] с разрешения правообладателя (Copyright

2017 AIP Publishing). б) Соответствующая картина температуры поверхности, полученная с помощью компьютерного

моделирования с использованием параметров из таблицы. в) График векторного поля скорости на поверхности капли,

полученный численно

Мы не рассчитывали инфракрасное изображе-

отличается от полученной численно в [25], где ячей-

ние с помощью уравнения (31) из [25], поскольку

ки БМ располагались только вблизи контактной ли-

считаем, что это уравнение неверно. В частности,

нии, и их количество увеличивалось с уменьшени-

уравнение написано в предположении, что тепловое

ем высоты капли. Зависимости расположения яче-

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

ек БМ от угла смачивания на рис. 1 и 5 похожи,

направлении. Кроме того, используется выражение

что естественно, поскольку соответствующие физи-

Стефана - Больцмана, при этом не учитывается, что

ческие параметры довольно близки.

максимум распределения Планка находится за пре-

делами диапазона длин волн, видимого инфракрас-

4. ЗАКЛЮЧЕНИЕ

ной камерой, т. е. инфракрасная камера пропускает

Мы выполнили трехмерное компьютерное моде-

значительную часть излучения.

лирование нестационарных неустойчивостей Маран-

На рис. 5 показаны вычисленные распределения

гони в испаряющихся лежащих каплях этанола и

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

силиконового масла капиллярного размера на на-

чений θ. Для большого контактного угла наблюда-

гретой подложке с высокой теплопроводностью. Мы

ется цветочная структура ячеек БМ. В то же вре-

проанализировали, насколько существенны различ-

мя ячейки занимают всю поверхность капли при до-

ные физические эффекты в этой задаче и обнару-

статочно малых углах смачивания. Такая ситуация

жили, что, помимо нестационарной трехмерной теп-

368

ЖЭТФ, том 159, вып. 2, 2021

Моделирование нестационарных неустойчивостей Бенара - Марангони. . .

Рис. 5. Структуры в распределениях температуры на поверхности капли, полученные при помощи компьютерного моде-

лирования капли этанола для θ = 29.2^◦ (а), 24^◦ (б), 20^◦ (в), 15^◦ (г)

лопроводности в капле, нестационарной трехмерной

ли, что уменьшение краевого угла приводит к увели-

динамики несжимаемой жидкости и диффузии па-

чению размера ячеек БМ, что в конечном итоге пре-

ров в воздухе, также необходимо учитывать поток

вращает цветочную пространственную структуру в

излучения между каплей и окружающей средой и

структуру, в которой ячейки БМ занимают всю кап-

теплопередачу в газовой фазе. Однако вклад таких

лю. Такое поведение наблюдалось эксперименталь-

явлений, как нестационарные эффекты диффузии

но в [27], но, насколько нам известно, до сих пор

паров в воздухе незначителен.

не было получено численно. Мы обнаружили, что

Интересно, что хорошо известное аналитическое

это поведение применимо и к испаряющейся капле

решение для скорости испарения в диффузионной

этанола [25]. Это устраняет несоответствие между

модели испарения [30], которое значительно упро-

поведением ячеек БM в [27] и в [25].

щает расчет теплопередачи в капле, также может

Наши результаты также согласуются с недавни-

быть использовано для расчета теплопроводности в

ми экспериментальными результатами [39], которые

газовой фазе.

показывают, что гидротермальные волны в капле

Наши результаты хорошо согласуются с экспе-

возникают только при относительно больших углах

риментальными результатами [25, 27]. Мы подтвер-

смачивания.

ждаем вывод этих работ о том, что неустойчивость в

данном случае представляет собой нестационарную

Мы полагаем, что полученные результаты пред-

конвекцию БМ, а не ГТВ.

ставляют собой полезный шаг на пути к лучшему

Также проясняется поведение неустойчивостей

пониманию термокапиллярной неустойчивости в ис-

БМ при уменьшении краевого угла. Мы обнаружи- паряющихся каплях.

369

ЖЭТФ, вып. 2

А. А. Гаврилина, Л. Ю. Бараш

ЖЭТФ, том 159, вып. 2, 2021

Финансирование. Работа выполнена при фи-

20.

M. K. Smith and S. H. Davis, J. Fluid Mech. 132,

нансовой поддержке Российского научного фонда

119 (1983).

(грант № 18-71-10061).

21.

M. K. Smith. Phys. Fluids 29, 3182 (1986).

22.

D. Schwabe, U. Möller, J. Schneider, and A. Schar-

ЛИТЕРАТУРА

mann, Phys. Fluids A 4, 2368 (1992).

D. Brutin and V. Starov, Chem. Soc. Rev. 47, 558

23.

R. J. Riley and G. P. Neitzel, J. Fluid Mech. 359,

(2018).

143 (1998).

R. G. Larson, AIChE J. 60, 1538 (2014).

24.

W.-Y. Shi, K.-Y. Tang, J.-N. Ma, Y.-W. Jia,

H.-M. Li, and L. Feng, Int. J. Therm. Sci. 117, 274

H. Y. Erbil, Adv. Colloid Interf. Sci. 170(1-2), 67

(2017).

(2012).

25.

S. Semenov, F. Carle, M. Medale, and D. Brutin,

D. Zang, S. Tarafdar, Yu. Yu. Tarasevich,

Appl. Phys. Lett. 111, 241602 (2017).

M. Du. Choudhury, and Tapati Dutta, Phys.

Rep. 804, 1 (2019).

26.

S. Semenov, F. Carle, M. Medale, and D. Brutin,

Phys. Rev. E 96, 063113 (2017).

F. Giorgiutti-Dauphiné and L. Pauchard, Eur. Phys.

J. E 41(3), 32 (2018).

27.

T.-S. Wang and W.-Y. Shi, Int. J. Heat Mass Transfer

131, 1270 (2019).

X. Shao, F. Duan, Y. Hou, and X. Zhong, Adv.

Colloid Interf. Sci. 275, 102049 (2020).

28.

L. Yu. Barash, T. P. Bigioni, V. M. Vinokur, and

L. N. Shchur, Phys. Rev. E 79, 046301 (2009).

R. Savino and S. Fico, Phys. Fluids 16, 3738 (2004).

29.

Н. Н. Лебедев, Специальные функции и их прило-

K. H. Kang, S. J. Lee, C. M. Lee, and I. S. Kang,

жения, Физматлит, Москва-Ленинград (1963).

Meas. Sci. Technol. 15, 1104 (2004).

30.

R. D. Deegan, O. Bakajin, T. F. Dupont, G. Huber,

H. Hu and R. G. Larson, J. Phys. Chem. B 110, 7090

S. R. Nagel, and T. A. Witten, Phys. Rev. E 62, 756

(2006).

(2000).

10.

W. D. Ristenpart, P. G. Kim, C. Domingues, J. Wan,

31.

H. Hu and R. G. Larson, J. Phys. Chem. B 106, 1334

and H. A. Stone, Phys. Rev. Lett. 99, 234502 (2007).

(2002).

11.

K. Sefiane, J. R. Moffat, O. K. Matar, and R. V. Cras-

32.

D. Ambrose and C. H. S. Sprake, J. Chem. Thermo-

ter, Appl. Phys. Lett. 93, 074103 (2008).

dyn. 2, 631 (1970).

12.

K. Sefiane, A. Steinchen, and R. Moffat, Colloids

33.

Н. А. Фукс, Испарение и рост капель в газообраз-

Surf. A 365, 95 (2010).

ной среде, Изд-во АН СССР, Москва (1958).

13.

D. Brutin, B. Sobac, F. Rigollet, and C. Le Niliot,

34.

Л. Д. Ландау, Е. М. Лифшиц, Теоретическая фи-

Exp. Therm. Fluid Sci. 35, 521 (2011).

зика, т. 6, Гидродинамика, Наука, Москва (1988).

14.

B. Sobac and D. Brutin, Phys. Fluids 24, 032103

35.

G. W. Kattawar and M. Eisner, Appl. Opt. 9, 2685

(2012).

(1970).

15.

F. Carle, B. Sobac, and D. Brutin, J. Fluid Mech.

36.

D. Thompson, L. Zhu, R. Mittapally, S. Sadat,

712, 614 (2012).

Z. Xing, P. McArdle, M. M. Qazilbash, P. Reddy, and

16.

K. Sefiane, Y. Fukatani, Y. Takata, and J. Kim,

E. Meyhofer, Nature 561, 216 (2018).

Langmuir 29(31), 9750 (2013).

37.

V. Fernández-Hurtado, A. I. Fernández-Dom´ınguez,

17.

G. Karapetsas, O. K. Matar, P. Valluri, and K. Se-

J. Feist, F. J. Garc´ıa-Vidal, and J. C. Cuevas, Phys.

fiane, Langmuir 28(31), 11433 (2012).

Rev. B 97, 045408 (2018).

18.

P. J. Sáenz, P. Valluri, K. Sefiane, G. Karapetsas, and

38.

S. V Patankar, Numerical Heat Transfer and Fluid

O. K. Matar, Phys. Fluids 26, 024114 (2014).

Flow, CRC Press (2018).

19.

J.-L. Zhu, W.-Y. Shi, and L. Feng, Int. J. Heat Mass

39.

T.-S. Wang and W.-Y. Shi, Int. J. Heat Mass Transfer

Transfer 134, 784 (2019).

148, 119138 (2020).

370