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

УДК 517.9

0 подходах к решению задачи о деформации межфазной границы в двухслойной системе с испарением *

В.Б. БекежановаО.Н. Гончарова2'3

1 Институт вычислительного моделирования СО РАН (Красноярск, Россия)

2 Институт теплофизики им. С.С. Кутателадзе СО РАН (Новосибирск, Россия)

3 Алтайский государственный университет (Барнаул, Россия)

On Approaches to Solving the Problem of an Interface Deformation in a Two-Layer System with Evaporation

V.B. Bekezhanova1, O.N. Goncharova2'3

1 Institute of Computational Modeling SB RAS (Krasnoyarsk, Russia)

2 Kutateladze Institute of Thermophysics SB RAS (Novosibirsk, Russia)

3 Altai State University (Barnaul, Russia)

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

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

БМ 10.14258/izvasu(2018)1-12

Mathematical modeling of the two-layer convective fluid flows is performed by the Oberbeck — Boussinesq equations and relations at the thermocapillary interface. Special attention should be paid to the problems with mass transfer at the interfaces. In particular, mass transfer can be a result of evaporation or condensation. Also, the problems with the vapor diffusion and thermodiffusion and diffusive thermal conductivity effects in the gas-liquid layers should be considered. The new solutions of the convection equations are presented to model two-layer flows. The flows are induced by the action of a longitudinal temperature gradient in the transversely directed gravity field. The interface assumed to be rectilinear when developing the analytical solution. The problem of finding the real interface position is solved with the help of the interface conditions. A derivation of the kinematic and dynamic conditions in terms of stream function and vorticity and an equation for the tangential velocity at interface is presented. A method of determination of a smooth deformable interface is discussed.

Key words: convective flows, evaporation, deformable surface, method of interface determination, exact solutions.

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

*Работа выполнена при финансовой поддержке РНФ (проект 15-19-20049)

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

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

Математическое моделирование конвективных течений включает построение новых точных решений определяющих уравнений и исследование влияния гравитационных и термокапиллярных эффектов, граничного теплового режима и способности границы абсорбировать испаряющуюся жидкость на характер течения с помощью полученных аналитических решений [2-4]. Система уравнений Обербека — Буссинеска, описывающая конвективные движения жидкостей, является достаточно сложной в силу своей нелинейности, высокого порядка и невозможности отнести ее к какому-либо классическому типу [5]. Построение решений дифференциальных уравнений специального вида позволяет перейти к задаче меньшей размерности. Точные решения дают возможность быстро и эффективно моделировать реальные течения жидкостей, в том числе течения в условиях фазового перехода на поверхностях раздела. Важность точных решений заключается также и в том, что их использование позволяет на качественном уровне выделить физические факторы, определяющие основные механизмы течения. В ситуации, когда задача является многопараметрической, именно точное решение позволяет изучить степень и характер влияния отдельных факторов и их взаимных комбинаций, уточнить математическую модель с целью адекватного описания исследуемых течений.

Моделирование двухслойных течений на основе точных решений в двумерном и трехмерном случаях, исследование устойчивости течений, которые могут быть описаны данными решениями, проведено в [2-4]. Новые точные решения, описывающие установившиеся конвективные течения с учетом испарения на термокапиллярной границе раздела в канале, получены в [4] для трехмерного случая при отсутствии осевой симметрии и в [1-3] — в рамках двумерной постановки. Указанные решения удовлетворяют уравнениям Навье — Стокса в приближении Обербека — Буссинеска, в которых дополнительно учитываются эффекты Соре (термодиффузии) и Дюфу-ра (диффузионной теплопроводности) [6-8] в газопаровой среде и которые являются аналогами решения Остроумова — Бириха [9-11].

При построении аналитического решения уравнений конвекции полагается, что граница

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

Постановка задачи и ее решение. Пусть система координат выбрана таким образом, что вектор силы тяжести g направлен противоположно оси Оу ^ = (0, —д)). Жидкость и газопаровая смесь заполняют бесконечные слои толщиной I и h. Ниже индексы г = 1, 2 будут определять верхний (газопаровой) г = 2 или нижний (жидкий) г = 1 слой системы. Представим двумерную задачу для определения компонент вектора скорости V = (и,ю), модифицированного давления р (отклонение давления от гидростатического значения), температуры Т и концентрации пара С, удовлетворяющих приближению Обербека — Буссинеска уравнений Навье — Стокса [8]:

У в |

0 X

-1

2Ю геометрия области течения

, П"

UiUix + ViUiy = —rjiPix + -¡От- Ащ

П"

+viviy = piy+Re Avi+(

Re nj Gr

R e2

(1)

Ti + ^C ),

Re2

(2) (3)

UiTix + ViTiy

+ Viy = 0,

-(ATi + Ö(AC)), (4)

nj

RePr-1

U2Cx + V2Cy = —(AC + a(AT )).

(5)

Здесь система уравнений представлена в безраз-

Пр Ga

мерной форме. При этом pт = рт +—i—ry. Кон-

i Re2

центрация пара в газопаровой фазе определяется уравнением диффузии (5). Отмеченные в уравнениях (2) и (4) слагаемые и уравнение (5) учитываются только для описания динамики и тепломас-сопереноса в верхнем слое системы. Здесь Re = u*h/vl — число Рейнольдса, Pr = v1/x1 — число Прандтля, Gr = ßlT*gh3/v2 — число Грасго-фа, Ga = gh3/v1 — число Галилея, g = |g|, Pe = u*h/D — число Пекле, nP = 1, Пр = 1/Р, П" = П" = v, nj = 1 nj = X, Р = P2/Pi, v = V2/V1, X = X2/X1, ß = ß2/ßi; u* — характерное значение скорости жидкости, h — характерный размер (например, высота слоя жидкости), Pi — относительное значение плотности, vt и Xi. — коэффициенты кинематической вязкости и температуропроводности, D — коэффициент диффузии пара в газе, ßi — коэффициент теплового расширения, y — концентрационный коэффициент плотности, коэффициенты Ö и а характеризуют эффекты Дюфура и Соре.

Точные стационарные решения системы уравнений (1) — (5) представимы в виде [1,2]:

4 3 2

ui = 24 L4 + \ L3 + Т С1 + Ус2 + с3,

о,

Ti = (ai + a2y)x + ßi, C = (bi + b2y)x + p y2

P\

2 d3 + у^ + dlj х + (6)

где дг и ф представляют собой полиномы седьмой степени, а ф — полиномы восьмой степени относительно поперечной координаты, коэффициенты Llj выражаются через физические параметры системы, сгк являются константами интегрирования, dlk выражаются через физические параметры, константы интегрирования с\ и коэффициенты агк ,Ьк.

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

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

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

Разработке численных методов для определения деформируемых поверхностей раздела, включая свободные границы, посвящены работы [13, 14]. В частности, в [14] построена математическая модель для исследования течений в системе тонких жидких слоев, находящихся под действием тепловой нагрузки. Идея построения модели основана на декомпозиции сложной задачи на унитарные элементы с некоторым набором правил их сопряжения друг с другом [15]. Каждый из модулей представляет собой однотипную локальную модель, в которой расчет неизвестных характеристик организован независимо от других модулей системы. При этом «набор правил», определяющий, например, общие границы между модулями и функциональные зависимости на них, соединяет эти модули в исходную систему. Определяющую роль при таком подходе играют физические и теплофизические свойства жидких слоев. Именно они управляют конвективными движениями в слоях и определяют положение и форму поверхностей раздела. Кроме того, модель можно расширять, добавляя к ней новые данные о физических процессах, обеспеченных выполнением законов сохранения. Характерной чертой разработанных методов является решение задачи в переменных «функция тока-вихрь», при котором отправной точкой расчетов оказывается решение уравне-

urv

u

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

Рассмотрим динамическое условие на границе раздела Г (размерный вид см. в [8]):

[Т • п] = 2аНп + Уга, (7)

где [Т • п] = Т11 г • п— Т21 г • п, Т — тензор напряжений Тг = —рг1 + 2рiviV(vi), V — тензор скоростей деформаций г-й жидкости, Уг = V — п(п • V) — поверхностный градиент, п — единичный вектор внешней (по отношению к 1-й жидкости) нормали, а — коэффициент поверхностного натяжения. В результате проектирования условия (7) на единичный вектор нормали к поверхности раздела (внешней по отношению к нижнему слою) п и единичный касательный вектор в получим в безразмерном виде условия:

2

—Р1 + п • D(vl)n— Не

— [—Р2 + НП ^ ^)П]= СНаК (8)

д©

2в • D(v1)n — 2р1?в • D(v2)n = —Ма-. (9)

д.

Предполагается, что имеет место линейная зависимость коэффициента поверхностного натяжения от температуры: а = 1 — МпСа(© — ©0). При формулировке условий использованы следующие обозначения: Не— число Рейнольдса, Са — капиллярное число (Са = ———1), Мп — число Ма-

рангони (Ма =

ат Т*

а0

© — безразмерная тем-

Р1 — Р2

0,

(10)

а кинематическое условие (12) удовлетворяется автоматически.

Если реальное положение границы раздела Г задается уравнением у = 7(х(в рассматриваемом случае 7 (х, ^ — это отклонение от положения у = 0), то единичный вектор нормали к границе раздела определяется следующим образом:

п = (п1, П2) = ( —

!х

1

лЛ+7Г г

(11)

Единичный касательный вектор будет определяться соответственно как в = ?2) = (п2, —п1), а кривизна кривой равна

К = 7хх/{1 + 72)3/2 (или К = 1/Н, Н — радиус кривизны в [13,17]). Зависимость функции 7 от «фиктивного» времени t при решении стационарной задачи позволяет организовать итерационный процесс расчета положения границы раздела.

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

VI • п = V2 • п = Vn, Уи

7г

л/г+72'

(12)

Здесь Уи — скорость движения границы раздела в направлении нормали п. Вихрь и функция тока связаны с компонентами вектора скорости следу-

ющим образом: ш = ух — —у, — = Фу

-Фх.

Касательная и нормальная скорости жидкости на границе раздела у, = V • в, уи = V • п выражаются через производные функции тока по направлению касательного и нормального векторов, соответственно:

дФ

дп

дФ

д.

(13)

пература (ее отклонение от относительного значения), К = 2Н — кривизна границы раздела Г (здесь плоской кривой); Н — средняя кривизна поверхности Г (относительно соглашения о знаке Н см. [8]), а0 и ат — константы, а0 — характерное значение поверхностного натяжения, ат — температурный коэффициент поверхностного натяжения (ат > 0 для большинства жидкостей).

Заметим, что при построении точного решения (6) условие (8) на границе раздела Г, задаваемой уравнением у = 0, выполняется в виде

Условия (8), (9), (12) могут быть записаны в терминах ш — Ф (см., например, [13, 16]). Но если соотношения (9), (12) записываются в терминах ш — Ф достаточно просто, то уравнение (8) предстоит продифференцировать по касательному направлению. Затем следует воспользоваться уравнениями (1), (2), записанными для обеих сред (здесь жидкости и газопаровой смеси), для вычисления производной функции давления по касательному направлению. Будем придерживаться подхода, разработанного в [13,17] и позволяющего переписать (8) в виде уравнения для касательной скорости у, на границе раздела, определяемой теперь уравнением у = 7(х,^).

Результат дифференцирования функций р1 и р2 по направлению касательного вектора в есть следствие применения уравнений (1), (2) и следующих соотношений

дш ду, дуи

(Дv,в) = ——, ((vV)v, в) = +уиди—шvи,

ду

дУи

(vt, в) = Ж — ^дт — уи7

д.

ду

п • D(v)n = —д? + УиК,

д'о

2в • D(v)n = 2—^ — ш + 2уаК. д.

(14)

(15)

(16)

у=

у

а

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

Таким образом, имеем:

дрх 1

= \-uit - (у1 ' ^)и1 + ЕАм1}в1 +

г , ._1 . Gт _ Ga -

+{-"и - (V1 • + ЖА"1 + Е2Т1 - Вё2}а2 =

Для функции 7, определяющей отклонение границы раздела от положения у = 0, условие (12) можно переписать в виде (см. (13)):

-у/1+72

0.

(20)

или

дФ

д"1- дv = --"1п

дЬ

да

- ю1пК7х\—

дщ- дvln

{"1: ^--+ Vln —--Ш1"1п} +

да дш

да От

Ga

+ - дт) + а2{ Ел Т - Е2 };

дп

Ее2

Ее2

(17)

др2 ^ РР

— = {-ри2ь - • \)и2 + ЕеАм2}а1 +

рР

+ {-РV2t - Р^2 • V)V2 + —-АV2 +

Ее

+Етв рТ2 - 1е2р+Ее трс Ь =

_ д"2: дv2n 2 к* Л

= -Р{ ~аГ~ - "2п^~ - "2пК7х}-

дЬ

_г дv2« д"2п - Р {"2: ~ + V2n-

да

да

да

рр , дш2 ) ш2"2п} + ЕеУ - дП2)

дп

, От — Оа Оа

+а2{ Ее22в рТ2- Е2р + Ее* 1рс}. (18)

В результате дифференцирования соотношения (8) вдоль границы Г (т.е. по касательному направлению) получим:

др1 др2 2 д , . . 2рр д . , , - И+д2+Ё-ед-а^^--кёТа (^)п)

1д

(Ка).

СаЕе да

Пусть в качестве характерной скорости выбрана скорость релаксации вязких напряжений. Тогда Ее = 1. Поскольку выполняется (15), то получим с учетом (17), (18) следующее уравнение для касательной скорости V::

(' - р> д+(1 -р)"- д- = 2(1 -рр д?+* • <19>

1 д к 1 / дш1 __дш2 \

Са да дп Р дп

-2(1 -РР ) д ("пК) + + (1 - Р)К7х"П + "п(ш1 - рШ2)-

-т2[(рв - 1)ОтТ + (1 - р)Оа + р^ОаС].

7г + ^1 + 7! = 0.

При выводе уравнения (19) мы воспользовались непрерывностью касательной скорости на границе раздела Г, а потому равенством на данной границе самих касательных скоростей ("1: = "2: = ":) и их производных по касательному направлению и по времени, а также равенством на границе раздела нормальных составляющих скорости ("1п = "2п = "п, см. (12)) и их производных по касательному направлению.

Динамическое условие (9) представимо в виде

шт - риш2 =

д( д" Мп— +2(1 - рр) ^ + 2(1 - рр)К": (21) да да

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

Общая схема решения задачи о нахождении поправки границы раздела. Заметим, что нормальная составляющая скорости "п, содержащаяся в правой части уравнения (19), вычисляется согласно (12): "п = Уп, а второе слагаемое в левой части этого уравнения можно заменить 1 д 2

на — — (":) . Для нахождения положения грани-2 да

цы раздела Г следует организовать итерационный процесс нахождения функции у = 7(Х^Ь), исходя из уравнений (19) и (20). Точное решение (6) продиктует необходимые для вычисления «правой части» V значения функций Т• С, а также = -щу (последнее потребуется для нахожде-дш

ния —— на границе раздела). дп

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

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

в [14].

п

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

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

Авторы выражают благодарность А.С. Овча-ровой за интерес к работе и плодотворные дискуссии.

Библиографический список

1. Гончарова О.Н., Резанова Е.В., Люлин Ю.В., Кабов О.А. Моделирование двухслойных течений жидкости и газа с учетом испарения // Теплофизика и аэромеханика. — 2015. — Т. 22, № 5.

2. Bekezhanova V.B., Goncharova O.N. Stability of the exact solutions describing the two-layer ows with evaporation at interface // Fluid Dynamics Research. — 2016. — Vol. 48, No 6 (061408).

3. Бекежанова В.Б., Гончарова О.Н., Резанова Е.В., Шефер И.А. Устойчивость двухслойных течений жидкости с испарением на границе раздела // Известия РАН. Механика жидкости и газа.— 2017. — № 2.

4. Bekezhanova V.B., Goncharova O.N. Study of the convective fluid flows with evaporation on the basis of the exact solutions in a three-dimensional infine channel // IOP Conf.Series Journal of Physics: Conference Series. — 2017. — Vol. 899, No. 032006. D0I:10.1088/1742-6596/899/3/032006.

5. Пухначев В.В. Нестационарные аналоги решения Бириха // Известия Алт. гос. ун-та. — 2011. — № 1/2 (69).

6. De Groot S.R. 1951 Thermodynamics of irreversible processes. — New York, 1951.

7. Gebhart B., Jaluria Y., Mahajan R.L., Sammakia B. Bouyancy-induced flows and transport — Berlin —Heidelberg —New York — London —Paris —Tokyo, 1988.

8. Andreev V.K., Gaponenko Yu.A., Goncharova O.N., Pukhnachov V.V. Mathematical

models of convection (de Gruyter Studies in Mathematical Physics). — Berlin/Boston: De Gruyter, 2012.

9. Бирих Р.В. О термокапиллярной конвекции в горизонтальном слое жидкости // ПМТФ. — 1966. — № 3.

10. Остроумов Г.А. Свободная конвекция в условиях внутренней задачи. — М. ; Л., 1952.

11. Шлиомис М.И., Якушин В.И. Конвекция в двухслойной бинарной системе с испарением // Гидродинамика. — 1972. — № 4.

12. Марчук И.В, Кабов О.А. Модель пленочной конденсации пара на криволинейных поверхностях // ДАН. — 2016. — Т. 466, № 1.

13. Овчарова А.С. Численное решение стационарной задачи Стефана в области со свободной границей // Вычислительные технологии.— 1999. — № 1.

14. Ovcharova A.S. Multilayer system of films heated from above // Int. J. Heat Mass Transfer. — 2017. — V. 114.

15. Воеводин А.Ф., Шугрин С.М. Методы решения одномерных эволюционных систем. — Новосибирск, 1993.

16. Пухначев В.В. Движение вязкой жидкости со свободными границами. — Новосибирск, 1989.

17. Ovcharova A., Stankus N. A Deformation and a Break of Hanging Thin Film under Microgravity Conditions // FDMP. — 2007. — Vol. 3,

No. 4.