МАТЕМАТИКА ИМЕХАН ИКА
УДК 517.9
Характеристики двухслойного течения с испарением в плоском канале принагревеснизу*
В.Б.впкежаноen 1, О.Н. Пончаунва2
1 Институт вычислительного моделирования СО РАН (Красноярск, Россия)
2 Ал тайский госсырственниайуниверситет(Барнтул,Россия)
Characteristics ofa Two-LayeaFlow with Evaporation in a Plane Channel Subjected to Heating from the Bottom
^Б. Beh^e^zhattova pO.N. GonchpwBa2
1 Insffiute ofComputational ModelmgSBRA0 (ргаероуыгsk.Russm) 2Ahai SteteUniversityiBarnaul, I^i^^^^aS
Исследуются особенностиконвективн!^ехр^:жип^ов, возюмнюнн« в деухслойнойсистемев усыивижы пвяеходс. Матимоымиеская модеть нля описвнияиснари-тальной 1^с^^лкцинрВвскон^^(^м^с^ризонт^йном ктнате овновааа нссспноотимымыоМТырОиоа—Буьан-еотоаивс-ноньй И ааво — (й^(^с^сти сонийонюанах нркермокапил-ляонвйфанынетазиклк.Тошгоорешекиеопреаелкн>щ1н-уОтвнениН — аншго г р ешднилОстроумова — Болиев — имеет групповую природу и позволяет учесть одновременное наличие горизонтального и вертикального гра-дяентовтемтсротярыи влипинетермоунффиеиоизтех э—фектси (пи^^й^с^с^о иобиасюиОяпарогиооввй смеси и на межфазной границе. Описаны тепловая и топологическая структуры совместных течений испаряющейся жидкости и смеси ее паров с инертным газом, изучены хснектеоисыикинеро сатзожтиря вверхидйтсон.Вфобо-те приведены новые результаты по исследованию устойчивости изучаемого точного решения уравнений конвекции. Представлены типичные формы возникающих характеристических возмущений в условиях равных продольных температурных градиентов на внешних стенках канала и ненулевого поперечного перепада температуры. Описаны определяющие механизмы, отвечающие за фор-мир оваттекаждоготипа структур.
Крнкеевгеескопк.конвектитныетечения, ттпаренит, тукнр1ерешетря, устойчивость, спектр характеристических возмущений.
БОН 0.14258/izvasu(2018)4-10
Введение. Ктмаыканамыы аычымня, сткнт-аеждиющився фиотаым кыныхтйтм «жидкость -пин», ш-нткн кннмымяюася а емыньыанкы, хи-
* Работа выполнена при финансовой поддержке РФФИ (проект 17-08-00291)
Features of convective regimes arising in atwo-layersystemwithaphase transition areinvestigeied. Amathemotical modeo tu Upsoe^tieOias eoeporctive oonvectioaév an inñnete herizrntcl chongelitiS^^sed on dhe hatear bred — Bou s sinepv coitus oxtmatton ohPhaNnvtec — StoVnrspuatéoht andon the relations an the thetmocapillasy ioSenogce. Anenest mlption ed the gooesningequatinps if thx COstpoumeg — Birikh solution analog. It has a group origination and allows one to take into account simultaneous presence ofhorizontalan— vrftiehltemyahatategfvUiryts med inflpance ofihetmodiOSusconaffects (disei:tand innerse) both in the gas-vapor mixture and on the interface. Thermal and topological patterns of the joint flows of an evaporating liquid and a mixture of its vapor with an inrrf s^are deteribsAChoeacieristict oeeopoe quahty in the upper layer are studied. New results on stability of the exact solution under given consideration are presented in the paper. Typical forms of arising characteristic perturbations are calculated for the case of e qual longitudin altemperatu re grad ients ontheexternal channelwalls andofnonzerottansversal tempetatoedrop. Goveemngmtehanisms responsible for the formation of each type of the structures are described.
Key words: Convective flows, evaporation, exact solutions, stability, spectrum of characteristic perturbations.
мической промышленности, материаловедении и биотехнологиях. Активное использование таких двухфазных систем требует детального изучения особенностей их функционирования, условий рационального применения и возможностей опти-
мизации. Одним из наиболее эффективных способов исследования разнообразных качеств жидкостных систем и характеристик возникающих конвективных режимов является метод математического моделирования, включающий построение корректных математических моделей, поиск, физическую интерпретацию и анализ свойств решений определяющих уравнений. Изучение характера влияния различных факторов на параметры тепломассопереноса позволяет предложить технические решения по повышению эффективности существующих жидкостных технологий и/или разработать принципиально новые подходы к решению задач термостабилизации в разнообразных системах с фазовым переходом, например таких, как микро- и мини-теплообменники, системы жизнеобеспечения, дистилляторы, двухфазные тепловые контуры и др.
Большинство аналитических работ, связанных с исследованием различных характеристик течений с испарением/конденсацией, выполняются в рамках математических моделей механики сплошной среды, а именно используют для описания испарительной конвекции уравнения На-вье — Стокса и переноса тепла и их аппроксимации. Дополнительной трудностью при таком подходе становится формулировка граничных условий на межфазной поверхности, которые должны учитывать процессы тепло- и массопереноса за счет испарения/конденсации через границу раздела сред [1,2]. Подробный обзор теоретических исследований по данной тематике представлен в [3]. Полезным и важным инструментом для изучения влияния различных тепловых, механических и физико-химических факторов на характер и интенсивность двухфазных течений являются точные решения уравнений конвекции. Такие решения позволяют на качественном уровне выделить основные физические механизмы, определяющие структуру основного течения, и исследовать характер и степень влияния отдельных факторов и их взаимных комбинаций. Особую ценность имеют решения, имеющие групповое происхождение. Именно решения групповой природы подразумевают сохранение свойств симметрии, заложенных при выводе основных уравнений и обеспечивающих тем самым физическое правдоподобие и реализуемость таких решений [4,5]. Уравнения На-вье — Стокса и Обербека — Буссинеска обладают богатыми групповыми свойствами [6], которые позволяют строить точные решения. При этом существуют классы решений, которые позволяют описать двухслойные течения с массопереносом через межфазную границу за счет испарения или конденсации в присутствии прямого и обратного термодиффузионных эффектов (эффектов Соре и Дюфура), как в плоской [7], так и в трехмерной геометрии [и]. Структура этих решений — обобще-
ний известного решения Остроумова — Бириха [9,10] — дает возможность проанализировать влияние различных факторов на характеристики режимов испарительной конвекции, включая свойства устойчивости, предложить классификацию типов течений и выделить типичные формы возмущений [11,12], а также апробировать разные типы граничных условий для функции концентрации пара в газовом слое (подробнее см. [3, 7,13] и список библиографии в них). Физическое правдоподобие и адекватность результатов, полученных на основе использования аналогов решения Остроумова — Бириха в задачах испарительной конвекции, подтверждается сравнениями с экспериментальными данными [14].
Характер и структура двухслойного течения и интенсивность эффектов испарения/конденсации определяются совместным действием разнородных факторов (естественная и термокапиллярная конвекция, гравитация, теплофизические свойства сред, управляющие воздействия, такие как расходы рабочих сред, тепловое воздействие на границах области течения и т.п., геометрия системы). Существенным образом возникающие режимы испарительной конвекции зависят от приложенной внешней тепловой нагрузки. В данной работе в рамках модели Обербека — Буссинеска исследуются характеристики совместного течения испаряющейся жидкости и парогазовой смеси в плоском горизонтальном канале, подогреваемом снизу, в условиях, когда на его стенках задан продольный градиент температуры. С помощью точного решения рассчитываются основные параметры режимов (скорости и температуры жидкости и парогазовой смеси, концентрация пара в газе, массовая скорость испарения), формы характеристических возмущений и нейтральные кривые. Проводится сравнение с решением аналогичной задачи, полученным в отсутствие поперечного перепада температуры.
1. Постановка задачи и вид точного решения. Пусть испаряющаяся жидкость и парогазовая смесь заполняют плоский горизонтальный канал с твердыми непроницаемыми стенками. Система координат выбрана так, что вектор массовых сил ё направлен противоположно оси Оу (ё = (0, -д)). Жидкий и парогазовый слой имеют толщины I и Ь соответственно. Межфазная граница раздела Г является термокапиллярной поверхностью, вдоль которой действуют касательные силы. При этом поверхностное натяжение а линейно зависит от температуры а = а0 — ж(Т — То); а0, Т0 — характерные значения поверхностного натяжения и температуры жидкости соответственно, ж > 0 — температурный коэффициент поверхностного натяжения.
Для описания совместного движения сред используется аппроксимация Обербека-Буссинеска
уравнений Навье — Стокса. Пар считается пассивной примесью, перенос пара в газе описывается уравнением диффузии. В газопаровом слое дополнительно принимаются во внимание эффекты Соре и Дюфура [15,16]. Таким образом, в стационарном случае определяющие уравнения имеют вид
(v •VW = -1 Vp' + vAv - g(/3T + yC), (1) P —
div v = 0, v •VT = x(AT + SAC), (2)
v •VC = D(AC + aAT), (3)
где подчеркнутые слагаемые и уравнение (3) учитываются только при моделировании течения в верхнем слое, v = (u,v) — вектор скорости, p' — модифицированное давление (отклонение p от гидростатического p' = p — р g • x, x = (x,y)), T — температура, C — концентрация пара в газе, р — относительное значение плотности, в — коэффициент теплового расширения, y — концентрационный коэффициент плотности, v и x — коэффициенты кинематической вязкости и температуропроводности, D — коэффициент диффузии пара в газе, коэффициенты S и а характеризуют эффекты Дюфура и Соре соответственно.
На внешних границах y = —l, y = h справедливы следующие условия:
Ul\y=-l = 0, U2\y=h = 0.
(4)
Т11 у=-1 = А1Х + Я-, Т2 | у=н = Л2х + Я+. (5)
Для концентрации пара задается условие нулевого полного потока пара
( дС дТ \
--к а-
V ду + ду )
0.
(6)
y=h
Считая, что поверхность Г остается недеформи-рованной в процессе движения (положим, у = 0), потребуем на ней выполнения кинематического и динамического условия. Первое выполняется тождественно, а второе может быть записано в проекциях на нормальный и касательный вектор как
(и1 (и2 дТ
Р1 = Р2, Р1и1~у = Р2^2~1--Ж дХ . (')
((у и у дХ
Условия непрерывности скорости и температуры на Г имеют вид
ul|y=0 = U2 |y=0,
Ч=о = T2U
(8)
Условие для тепловых потоков и соотношение баланса массы пара, включающие влияние термодиффузии и диффузионной теплопроводности на перенос тепла и массы, а также уравнение для концентрации насыщенного пара на границе раздела имеют вид
дТ1 дТ2 е dC
--К2—--SK2 —
ду ду ду
-LM,
(9)
ГЛ (дО" дТ M = —Dp2[— + а— v ду ду
= C 11 + ЛТ^
y=0
си = с* [1+ЧЧ=о -4. (10)
Здесь к — коэффициенты теплопроводности, Ь — скрытая теплота испарения, М — массовая скорость испарения, е = Ьу/^*Т2), у — молярная масса испарившейся жидкости, R* — универсальная газовая постоянная, С* — концентрация насыщенного пара при Т2 = Т0. Здесь и всюду ниже индексы j = 1 и j = 2 используются для идентификации характеристик среды в жидком и парогазовом слое соответственно. Для замыкания задачи дополнительно задается расход парогазовой смеси
/12
R =/Р2'и2 (у) (у. (11)
о
Задача (1) — (11) допускает аналог решения Остроумова — Бириха [17]
= uj (У), vj = 0 P'j = P'j (x,У),
Tj = (ai + а2у)х + (y),
C = (bi + Ь2у)х + ф(у).
(12)
Подробная постановка задачи с обсуждением типа граничных условий для температуры и концентрации пара на твердых границах, ее обоснование и анализ, описание алгоритма нахождения неизвестных функций и констант интегрирования даны в [3,13,14].
В силу второго условия из (8) справедливо равенство а1 = о\ = А, и распределение температуры в слоях имеет вид
Tj = (A + а2у)х + fy(у), j = 1, 2.
(13)
y=0
Решение (12) допускает случай, когда продольные градиенты температуры на стенках и границе раздела равны, т. е. А1 = А2 = А [13]. При этом, если Я- = Я+, то на границах канала прикладывается равная тепловая нагрузка. Условие Я- > Я+ соответствует ситуации, когда нижняя стенка канала горячее верхней (равномерный нагрев снизу), при Я- < Я+ — нижняя стенка холоднее верхней (равномерный нагрев сверху). Равенство градиентов А1 и А2 обеспечивает применимость метода нормальных мод для исследования линейной устойчивости рассматриваемого точного решения, который и будет использован ниже. В настоящей статье рассматривается случай Я- > Я+ и изучается влияние поперечного перепада температур на критические параметры устойчивости и типы наиболее опасных возмущений.
2. Задача об устойчивости точного решения. Исследуем линейную устойчивость точного решения (12) относительно малых нестационарных возмущений.
u
j
2.1. Безразмерные параметры. Выберем в качестве характерных масштабов длины, скорости, давления и температуры следующие величины: Ъ, V2/h, р^у^/Ъ?, д- соответственно. Тогда безразмерные переменные пространства суть £ = х/Ъ, п = у/Ъ, а времени — т = v2t/ЪЪ2. Для каждого параметра среды ш- определим безразмерный аналог ш- = ш-/ш2. Тогда индексу 2 = 1 соответствует область —1/Ъ < п < 0, индексу 2 = 2 — область 0 < п < 1.
При выбранном способе обезразмеривания исходная задача характеризуется следующими безразмерными параметрами и критериями подобия: Gr = дв2д-Ъ3/^^, Рг = v2/x2, Ga = дЪ/2, Le = D/x2, О = АЪ2/д-. Здесь Сг, Рг, Са, Ье — числа Грасгофа, Прандтля, Галилея и Льюиса, Q — параметр тепловой нагрузки.
2.2. Спектральная задача. Рассмотрим малые двумерные нестационарные возмущения скорости и- (£,п,т) = (и (£,п,т ),У-(£,П,т)), давления Р'(£,п,т), температуры ©'(£,п,т) и концентрации S'(£,п,т) решения (12). Предположим, что эти функции пропорциональны ехр [г (ах£ — Хт)] и определяют нормальные волны с комплексным декрементом Х = Хг + гХ4, описывающим эволюцию возмущений во времени, и волновым числом ах вдоль оси £.
Линеаризация уравнений (1)—(3) вблизи стационарного решения (12) приводит к системе для амплитуд малых возмущений («штрихи» у безразмерных параметров и функций опущены):
—1/Ъ <п< 0 : гахи1 + V/ = 0,
—гхи/ + гахи1и1 + и'/у1 =
= — ^ р/ + v(U'^ — аХи0 , Р
—гХУ/ + гаХи1У1 =
= — 1 Р/ + V (V/' — а2ХУ/) + ^г©', Р
—гХ©/ + гахи/©/ + иТ^ + VlTlv =
= РГ (©'/' — а2©/) ,
0 <п< 1: гахи2 + V' = 0, —гХи2 + гах^Щ + ^У = = —гахР2 + и2 — а.х и2, —гХУ2 + га.хи2У2 = = —Р2 + У2 — а2хУ2 + Gr©2 + 7 GaS, —гХ©2 + га.хи2©2 + и2Т2£ + У2Т2П =
©2' — ах©2 + 6 (С'' — ах С)
(14)
(15)
1
РГ
—iХS + гаxU2S + + У2Сп =
= рг(&' — а2xS + ад-(©'' — а1 ©)).
На твердых стенках и поверхности раздела для амплитуд возмущений заданы следующие граничные условия:
П = —1/Ъ : и' = У/ = ©' = 0,
п
1 : и2 = У2 = ©2
S' + ад+ ©2 = 0,
п= ©/
0 :
= ©2,
и' = и2, Р/ — Р2
У' = У2
0,
2 ^рУ/ — У2),
Ма
(16)
и2 — vрU'1 + гах(У2 — У) = "о" гах©,
к©/ — ©2 — = ^^ + ад-©').
2 д- ^д- у 7
Здесь © определяет общую температуру обеих сред на Г (© = ©-, 2 = 1,2), Ма = жАЪ2/^2р2) — число Марангони. Теперь «штрихи» в уравнениях (14)-(16) обозначают дифференцирование по переменной п. При выводе условий (16) предполагалось, что возмущения функций (12) не приводят к возмущению межфазной границы, т.е. Г остается недеформированной.
Итак, соотношения (14) - (16) составляют спектральную задачу для комплексного декремента Х, при этом амплитуды возмущений суть неизвестные функции, определяющие характеристические возмущения типа нормальных волн. Полученная задача решалась численно методом ортогонали-зации [18], который был адаптирован для случая системы с внутренней границей раздела [19]).
2.3. Карта режимов неустойчивости. Механизмы неустойчивости. Основная проблема состоит в определении критических тепловых нагрузок, приложенных на внешних границах канала, и типа наиболее опасных возмущений, которые приводят к потере устойчивости.
На рис. 1 представлена карта режимов неустойчивости в плоскости (ах, Ма), полученная для пары сред НРЕ-7100-азот со следующими физическими параметрами ш = {ш/,ш2}: V = {0.38• 10-6,0.15• 10-4} м2/с, р = {1.5• 103,1.2} кг/м3, в = {1.8 • 10-3,3.67 • 10-3} К, к = {0.07, 0.02717} Вт/(м-К), х = {0.4 • 10-7,0.3 • 10-4} м2/с, ж = 1.14 • 10-4 Н/(м^К), Ь = 1.11 • 105 (Вт-с)/кг, ц = 0.25 кг/моль, D = 0.7 • 10-5 м2/с, 7 = —0.5, С, = 0.45, 6 = 10-5 К, а =5 • 10-3 К- . Толщины нижнего и верхнего слоя принимались равными I = 3 и Ъ =5 мм соответственно, расход R = 9.6 • 10-6 кг/(м-с), относительная температура То = 20°С, д = 9.81 м/с2, д- = 20°С. Области неустойчивости лежат внутри нейтральных кривых 1 и 2, построенных для случаев нулевого (д+ = 20°С) и ненулевого (д+ = 17°С, подогрев снизу) перепада температур соответственно. Знак величины Ма совпадает со знаком продольного градиента А, так, что при Ма > 0 стенки канала нагреваются в направлении оси Ох,
при Ма < 0 — охлаждаются. На карте в области соответствующих волновых чисел указаны типы наиболее опасных возмущений, приводящих к потере устойчивости основного течения: К — конвективные ячейки, С — структуры смешанного типа, В — вихревые структуры, Т — термокапиллярные структуры. Индексы «+» и «— » указывают на режим тепловой нагрузки, при котором формируются структуры соответствующего типа: при нагреве (Ма > 0) или охлаждении (Ма < 0) стенок.
Характерной особенностью течений в канале с равной тепловой нагрузкой является термокапиллярная неустойчивость, при которой формируются только возмущения термокапиллярного типа Т1 и Т2 (рис. 1). При потере устойчивости, сопровождающейся образованием структур Т1, в системе возникают чередующиеся тепловые пятна (темные области соответствуют холодным зонам, светлые — горячим) с ядрами на границе раздела. Коротковолновые возмущения приводят к образованию термокапиллярных структур Т2 с шахматной упаковкой тепловых пятен, верхний ряд которых локализован на Г, нижний — внутри жидкого слоя. Заметим, что основное течение при равной тепловой нагрузке является течением термокапиллярного типа с возвратным движением в жидком слое (типичные распределения скорости и температуры показаны на рис. 2(а), подробнее классификацию течений см. в [7]).
Нагрев снизу имеет дестабилизирующий эффект, обусловленный неустойчивой температурной стратификацией в системе (распределения скорости и температуры при нагреве снизу показаны на рис. 2(б)), значительно расширяет область неустойчивости и приводит к смене типов возникающих характеристических возмущений. Потеря устойчивости под действием длинноволновых возмущений сопровождается появлением конвективных ячеек, при этом форма гидродинамических возмущений не зависит от типа тепловой нагрузки (нагрев или охлаждение
Ма
-52.5 -105
0 2 4 6 8 а,
Рис. 1. Карта режимов неустойчивости: нейтральные кривые Ма(ах), 1 — д- = д+ = 20оС, 2 — д- = 20оС, д+ = 17°С; типы и формы наиболее опасных возмущений
стенок, определяется знаком Ма) и представляет собой вихри с попарно противоположной циркуляцией. Форма тепловых возмущений может меняться в зависимости от знака Ма. При Ма > 0 возникают чередующиеся тепловые пятна, деформированные основным течением (структуры К +, рис. 1). При Ма < 0 тепловые пятна локализованы в приповерхностной зоне жидкого слоя и за счет интенсивной конвекции переносятся поперек парогазового слоя (структуры К-, рис. 1).
С увеличением ах конвективные ячейки сменяются структурами смешанного типа (структуры С, рис. 1), которые отвечают сосуществованию конвективного и термокапиллярного механизмов неустойчивости. Тепловые пятна в жидкости расщепляются и имеют по два ядра. Верхние ядра сносятся основным течением за счет термокапиллярного растекания, а нижние дрейфуют в центральной части жидкого слоя и возникают благодаря конвективному механизму, поскольку основное течение характеризуется формированием горячего термоклина внутри жидкости (рис. 2(в)).
Коротковолновые возмущения приводят к воз-никовению вихрей с попарно противоположной циркуляцией, при этом ядра вихрей совмещены с ядрами тепловых возмущений (структуры В, рис. 1). Подобный тип неустойчивости реализуется при внешней тепловой нагрузке, обеспечивающей устойчивую температурную стратификацию в жидком слое (рис. 2(г)), которая, однако, не гарантирует устойчивость основного течения. Имеет место термокапиллярное движение жидкости от горячего полюса к холодному, которое в силу неразрывности среды порождает движение в вертикальном направлении, а следовательно, к формированию вихрей.
О х, мы 50 О х, мы 50
Рис. 2. Распределение скорости (сплошная линия) и температуры в двухслойной системе: (а) — А = 3°С/м; (б) — А = -4.48°С/м; (в) — А = 7.16°С/м; (г) — А = -9.82°С/м
Заключение. С помощью точного решения задачи испарительной конвекции исследовано влияние поперечного перепада температур на характеристики устойчивости совместного течения тонкого слоя жидкости и парогазовой смеси в условиях фазового перехода в плоском горизонтальном канале. Рассчитаны критические тепловые нагрузки, при которых происходит потеря устойчивости. Даже незначительный нагрев снизу существенно дестабилизирует основное течение, меняя основной механизм неустойчивости. При равной тепловой нагрузке на внешних границах канала определяющим фактором является эффект Марангони, приводящий к формированию термокапиллярных структур разных типов в жидком слое. В условиях ненулевого поперечно-
го перепада проявляется конвективный механизм неустойчивости, вызывающий появление разнообразных структур. Чисто конвективная неустойчивость проявляется возникновением конвективных ячеек с разной локализацией тепловых возмущений в зависимости от характера тепловой нагрузки на стенках (нагрев или охлаждение в направлении продольной оси). В условиях сосуществования термокапиллярного и конвективного механизмов наблюдаются структуры смешанного типа. При тепловых нагрузках, обеспечивающих устойчивую температурную стратификацию в жидком слое, возникают вихревые структуры, индуцированные термокапиллярным движением жидкости на поверхности раздела из горячей области в холодную.
Библиографический список
1. Margerit J. et al. Interfacial nonequilibrium and Benard-Marangoni instability of a liquid-vapor system // Phys. Rev. E. - 2003. - T 68 (041601).
2. Das K.S., Ward C.A. Surface thermal capacity and its effects on the boundary conditions at fluid-fluid interfaces // Phys. Rev. E. — 2007. -T 75 (065303).
3. Бекежанова В.Б., Гончарова О.Н. Задачи испарительной конвекции (обзор) // Прикладная математика и механика. — 2018. — Т. 82. — T 2.
4. Андреев В.К., Капцов О.В., Пухна-чёв В.В., Родионов А.А. Применение теоретико-групповых методов в гидродинамике. — Новосибирск, 1994.
5. Андреев В.К., Бублик В.В., Бытев В.О. Симметрии неклассических моделей гидродинамики. — Новосибирск, 2003.
6. Пухначеев В.В. Симметрии в уравнениях Навье-Стокса // Успехи механики. — 2006. — T 6.
7. Bekezhanova V.B., Goncharova O.N., Shefer I.A. Analysis of an exact solution of problem of the evaporative convection (Review). Part I. Plane case //J. Sib. Fed. Univ. Math.&Phys. — 2018. — Vol. 11, T 2.
8. Bekezhanova V.B., Goncharova O.N., Shefer I.A. Analysis of an exact solution of problem of the evaporative convection (Review). Part II. Three-dimensional flows // J. Sib. Fed. Univ. Math.&Phys. — 2018. — Vol. 11, T 3.
9. Бирих Р.В. О термокапиллярной конвекции в горизонтальном слое жидкости // ПМТФ. — 1966. — T 3.
10. Остроумов Г.А. Свободная конвекция в условиях внутренней задачи. — М. ; Л., 1952.
11. Родионова А.В., Резанова Е.В. Исследование устойчивости двухслойного течения жидкости // ПМТФ. — 2016. — T 4.
12. Резанова Е.В., Шефер И.А. О влиянии тепловой нагрузки на характеристики течения с испарением // Сиб. журнал индустриальной математики. — 2017. — T 2 (70).
13. 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, T 6 (061408).
14. Гончарова О.Н., Резанова Е.В., Люлин Ю.В., Кабов О.А. Моделирование двухслойных течений жидкости и газа с учетом испарения // Теплофизика и аэромеханика. — 2015. — Т. 22, T 5.
15. De Groot S.R. Thermodynamics of irreversible processes. — New York, 1951.
16. Gebhart B., Jaluria Y., Mahajan R.L., Sammakia B. Bouyancy-induced flows and transport. — Berlin — Heidelberg — New York — London — Paris — Tokyo, 1988.
17. Шлиомис М.И., Якушин В.И. Конвекция в двухслойной бинарной системе с испарением // Гидродинамика. — 1972. — T 4.
18. Годунов С.К. О численном решении краевых задач для систем линейных обыкновенных уравнений // УМН. — 1961. — Т. 16, вып. 3 (99).
19. Гапонов С.А., Маслов А.А. Развитие возмущений в сжимаемых потоках. — Новосибирск, 1980.