О влиянии зависимости вязкости от температуры на стационарные конвективные течения в ячейке Хеле-Шоу Текст научной статьи по специальности «Физика»

Механика

УДК 532.5.032 DOI: 10.14529/mmph170206

О ВЛИЯНИИ ЗАВИСИМОСТИ ВЯЗКОСТИ ОТ ТЕМПЕРАТУРЫ НА СТАЦИОНАРНЫЕ КОНВЕКТИВНЫЕ ТЕЧЕНИЯ В ЯЧЕЙКЕ ХЕЛЕ-ШОУ

В.А. Демин, М.И. Петухов

Пермский государственный национальный исследовательский университет, г. Пермь, Российская Федерация E-mail: demin@psu.ru

Представлены результаты численного моделирования стационарных конвективных течений в ячейке Хеле-Шоу при подогреве снизу. Выполнен линейный анализ устойчивости механического равновесия и проанализирован переход от одновихревого течения к двухвихревому при учете в математической модели зависимости вязкости жидкости от температуры. Для различных вариантов тепловых граничных условий получены поля функции тока при разных числах Рэлея, по которым определены сценарии интенсификации конвекции с ростом надкритичности.

Ключевые слова: ячейка Хеле-Шоу; подогрев снизу; стационарные режимы; неоднородность вязкости.

Введение. В работе проведено численное исследование влияния неоднородности вязкости и тепловых граничных условий на стационарные конвективные движения в вертикально ориентированной ячейке Хеле-Шоу при равномерном подогреве снизу. В рассматриваемых условиях при небольших перепадах температуры в жидкости имеет место механическое равновесие, а при малых надкритичностях возникают стационарные одно- или двухвихревое конвективные движения жидкости, область существования которых и форма исследовались в данной работе. Впервые обстоятельное экспериментальное наблюдение различных ламинарных конвективных режимов в подогреваемой снизу вертикальной ячейке Хеле-Шоу проводилось в работе [1]. Линейный анализ устойчивости течений в случае идеально теплопроводных широких граней был выполнен в [2]. Однако, из-за специфических тепловых условий на широких гранях и выбранного соотношения сторон 2:10:20, последовательность регулярных режимов оказалась достаточно простой. После прохождения порога устойчивости равновесия в системе сначала рождался одновихревой стационарный режим. По мере увеличения надкритичности ось симметрии основного вихря начинала отклоняться, а скорость вращения расти. В определенный момент система проходила точку бифуркации, в которой рождался автоколебательный четырехвихревой режим с перезамыканием угловых вихрей.

Позднее последовала серия экспериментальных и теоретических работ [3-5], в которых использовалась полость с соотношением сторон 2:20:40 и было обнаружено значительно большее разнообразие как переходных, так и устойчивых конвективных течений. В частности был обнаружен и описан новый, так называемый пульсационный режим, колебательное поведение которого проистекает из неустойчивости стационарного двухвихревого течения, когда спонтанно один вихрь начинает доминировать над другим. В какой-то момент большой вихрь целиком поглощает маленький, и течение на время становится одновихревым. Затем на фоне одновихревого течения в одном из нижних углов полости рождается небольшой вихрь противоположной закрутки. Этот вихрь со временем растет, достигает в размерах примерно 1/3 высоты полости, и, взаимодействуя с основным течением, в некоторый момент времени опять поглощается им. Через некоторый промежуток времени этот процесс повторяется. В итоге наблюдаются пульсации нижнего вихря на фоне основного одновихревого течения.

Эксперименты, проведенные на трансформаторном масле, привели к нетривиальному результату. Оказалось, что автоколебательное четырехвихревое течение характеризуется разными по величине верхними и нижними вихрями, что было объяснено в [5] наличием у данной жидкости сильной зависимости вязкости от температуры. В тоже время в отмеченных работах не обра-

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

d £ z

l

2,4 Ra

1,6

0,8

0

у / /

■ ' , - . '

12

31 22

21 11

0 0,1 0,2 а 0,3

а) б)

Рис. 1. Ячейка Хеле-Шоу (а). Пороги устойчивости механического равновесия для первых пяти гидродинамических мод (б) для двух значений е. Сплошные и штриховые линии отвечают, соответственно, однородной и неоднородной вязкости

1. Постановка задачи. Ячейка Хеле-Шоу представляет собой полость в форме прямоугольного параллелепипеда, одна из сторон которого много меньше двух других (рис. 1, а). По причинам, упомянутым выше, будем рассматривать ячейку с соотношением сторон 2С:1:к = 2:20:40; С,I, к - размерные толщина, длина и высота полости. Изучим конвективное поведение рабочей жидкости, находящейся в условиях равномерного нагрева снизу при наличии поля тяжести. При этом предполагается, что вязкость уменьшается с ростом температуры [6] по линейному закону V = v0(1 - еТ), где е- размерный коэффициент зависимости вязкости от температуры.

В качестве граничных условий рассмотрим два случая: 1) теплоизолированные и 2) идеально теплопроводные узкие вертикальные грани. На верхнем и нижнем теплообменниках поддерживалась определенная разность температур. Широкие грани обладают конечной теплопроводностью, так чтобы имелось максимальное соответствие эксперименту, в котором полость была ограничена плексигласовыми пластинами толщиной 2 ■ 3 см.

Для численного моделирования конвективных течений использовалась система уравнений тепловой конвекции в приближении Буссинеска [6] при наличии неоднородности вязкости:

^ + ^= -Ур + (1 -кТ)Дv - 2к(УТУ^ -кУТхюХь + ЯаТу, (1)

7)Т

Pr— + (vV)T = AT, divv = 0 . dt

(2)

Здесь V, р, Т - безразмерные поля скорости, давления и температуры, у - единичный вектор, направленный вертикально вверх. Система уравнений содержит безразмерные управляющие параметры: числа Рэлея и Прандтля, а также параметр, характеризующий зависимость вязкости от

температуры: Яа = gp0d3/v0c, Рг = v0/с, к = е0, где v0 - коэффициент кинематической вязкости в реперной точке; х, в - коэффициенты температуропроводности и теплового расширения, 0 - характерная разность температур на расстоянии, равном полутолщине слоя С, g - величина ускорения свободного падения. Дополнительными безразмерными параметрами являются длина и

h

x

высота полости: L = l/d , H = h/d . За единицы измерения расстояния, времени, скорости, давления и температуры выбраны, соответственно: d, d2/v0, j/d, pv0j/d2, 0 (p - средняя плотность жидкости). Число Прандтля в расчетах для простоты принималось равным Pr = 6, что приближенно соответствует воде.

Система уравнений (1)-(2) предполагает использование приближения плоских траекторий, следовательно, можно ввести функцию тока vх = Э¥/Эу , v = —Э¥/Эх. После исключения давления и получения уравнений в терминах «вихрь - функция тока» решение системы искалось в виде ¥ = у (x,y,t) cos (pz/2). На границах полости выполнялось условие прилипания жидкости. У поля температуры выделялась равновесная линейная часть: T = — y + T', где T' = в(x,y,t)(1 + 2ap_1cos(pz/2)). Граничное условие для температуры на широких гранях имеет

вид: ЭТyЭz|г =— a(T' — T0), где a = XJXp - безразмерный коэффициент теплоотдачи, 1m, If -

коэффициенты теплопроводности массива и жидкости, соответственно, Р - безразмерная толщина массива, T0 - температура окружающей среды, которая в данной задаче бралась за начало отсчета. Уравнения усреднялись по z, после чего подлежали определению у и 9, т.е. поля амплитуд функции тока и отклонения температуры от линейного профиля. Уравнения для амплитуд и начально-краевые условия имеют вид:

8 (ЭуЭф ЭyЭjN Эt 3pPr V Эх Эу Эу Эх

f

+

1 + ky — k I 1 +

16a 3P

Л

в

Э2у Э j p

2

Эх2

+

2—Т j

+k

V

1+-

16a

p

((Э2в Эв

+

V

Эх2 Эу2

j+-

p

2

2 —a 3

3p

Эх Эх Эу Эу ,

Эу Эв + Эу Эв

Эх Эх Эу Эу

16a Л(

Эу p2 Эу

Т "эУ

+

■2kI1+

2

3p2

djdq+djdq

Эх Эх Эу Эу

2п Л

Э j + Э2у Э2в Э уЭ2в Э2у Э2в

3p2 + 16a Эу ' Эх2 Эх2 ' Эу2 Эу2

' 8a2 Л

1 + a +—

Pr I 1 +

2a(4+a)

p

Эв

+

Эt

1+

V

3p

ЭхЭу ЭхЭу Эу Эв ЭуЭвЛ

^ Ra Эв, p Эх

(3)

2a (4 + a) Л( Э2в Э2вЛ

p

+

Эх Эу Эу Эх a(2+a)

2 + a Эу p Эх

+

/ v

>2,„ -ч2

Эх2 Эу2

2

в,

(4)

Эу

(5)

(6)

Э2у Э 2у

где введено следующее обозначение р = —— +--— ;

Эх Эу

г = 0: у= р= в= 0;

у = 0, Н: в = 0, у = У = 0; х = 0, Ь: у = ^ = 0;

Эх Эу

идеально теплопроводные узкие вертикальные грани: в = 0;

эв п

теплоизолированные узкие вертикальные грани: — = 0 .

Эу

2. Линейный анализ устойчивости. Чтобы оценить влияние коэффициента теплоотдачи на течение, а также для выбора наиболее близкого к эксперименту [4] значения этого параметра (сравнение велось по пороговому значению моды 1,1), сначала была аналитически решена задача линейной устойчивости механического равновесия. Для произвольных к граница устойчивости имеет вид:

Ra nm =

_4 т2

p L

, H. 1 +—k 2

2

n Л ( m

I) + (m

V

2 Л2 (/ \2

n Л ( m

I) + [ H

J V

1

+ — 4

л

X

4

2

2

n

X

(

p + 4a

2 (2 + a)

+

2+a

i+1

4

n )2 ( m

l ) +1 m

2 )

-1 ))

J J

Здесь т и п - целые числа, определяющие периодичность возмущений по осям х и у. Найденные значения критического числа Рэлея в пределе к = 0 совпадают с известными случаями теплоизолированных (а = 0) и идеально теплопроводных широких граней (а = ¥). Графически значения первых пяти гидродинамических мод в исследованном диапазоне коэффициента теплоотдачи представлены на рис. 16. Все дальнейшие расчеты проводились при а = 0,2, что примерно соответствует плексигласу при толщине стенок 2 см. Пороговым значением числа Рэлея для данного коэффициента теплоотдачи является Яа ~ 0,9.

3. Метод решения нелинейной задачи. Для анализа полных нелинейных уравнений (3)-(4) применялся метод конечных разностей. Краевая задача (3)-(6) в терминах конечных разностей с учетом обозначений школы А.А. Самарского [7, 8] записывается в форме

j

8

3pPr

(

У о jo

V x y

У o jo

y x J

+

1 + ky - k I 1 +

16a 3P2

J

q

i

+k

1+

16a

(qxx +qyy )j +

2

jxX + jyy ■ ))

2)

-j

P

+

2

—a 3

( )) jd + Уово +Уово

x X y y Jj

Уодо + y о q

x x y y J J

(

p 4

У о

■2kI 1 +

16a

jodo +jodo

yy

3p

2

\

3p + 16a

jo + Уд + Уд + 2yo о до

2 (2 + a)

Ra до

(

Pr

1+

2a (4 + a)

p

qt =-

г

1+a+

3p

x y x y J

2 )(

У о до - У о до

V x y y x

p

(7)

2+a

p

-У +

+

1+^ 1 д+dy)-^ q,

p J 2

(8)

í,N 2-11

r¡n _ r¡n N1, j ~ N1-1, j ■

(9) (10)

у = 0, H: = 0, eiN2 = 0, yn y = 0, jji = lyl/hl, jlN2 = 2yi

x = 0, L: yj = yNi,j = 0, jj = 2yn, , j} = 2yN1_1,} /hx2 ;

теплопроводные узкие вертикальные грани вЩ/ = 0, вгЩ1 j = 0 ;

теплоизолированные узкие вертикальные грани вЩ j = вЩ /

Здесь j = y - + yyy - с точностью до знака вихрь скорости; i, j - номера узлов, а N1, N2 - их

число по x и у. Для jrn границах использовались формулы Тома [8]. Введенные обозначения для конечных разностей отвечают следующим односторонним и центральным разностям

rn+1 _ rn rn _ rn rn _ rn

r _ J i, j J i, j r _ J i+1, j J i-1, j r _ J i+1, j J i, j r _

Jt , J o ~ , , J x , , J x

nn J i, j J i-

1, j

Т х 2кх кх кх

Здесь Т, кх, ку - шаги по времени и координатам х, у. Таким образом, алгоритм был разработан в соответствии с явной схемой решения уравнений в частных производных и основан на двухполе-вой методике [8]. При аппроксимации производных по времени и производных по координатам использовались, соответственно, односторонние и центральные разности. Расчет начинался с числа Рэлея, соответствующего равновесию. Таким образом, система стартовала из состояния покоя. В ходе пошагового изменения числа Рэлея в качестве начальных условий использовались поля в предыдущий момент времени, т.е. применялся метод продолжения по параметру. Это позволяло имитировать эксперимент, в котором постепенно после фиксации установившегося режима последовательно малыми шагами увеличивалась разность температур на теплообменниках.

a

y

xx

y

x

x

Шаг по времени выбирался из соображений устойчивости численной процедуры. При решении уравнения Пуассона для функции тока в каждый момент времени применялся метод простых итераций. Рабочее количество узлов в сечении канала было равно 25x49. В ходе расчетов использовался метод установления. Компьютерный код был реализован на языке программирования Б0ЯТЯЛК-90. Система уравнений (7), (8) с граничными условиями (9), (10) решалась на суперкомпьютере «ПГУ-Тесла» Научно-образовательного центра Пермского государственного национального исследовательского университета «Параллельные и распределенные вычисления».

4. Результаты расчета в случае однородной вязкости. Сначала проанализируем результаты численного моделирования стационарных течений в ячейке Хеле-Шоу в пределе е = 0 (однородная вязкость). Пусть для определенности узкие вертикальные грани теплоизолированы.

45-

5 10 15 20 25 5 10 15 20 25 5 10 15 20 25

Рис. 2. Поля функций тока стационарных конвективных течений для случая однородной вязкости, а) -однових-ревой режим при Ра = 1,6; б) и в) - соответственно двухвихревые режимы с подъемным и опускным течением с

середине полости при Ра = 2,7

В соответствии с линейным анализом устойчивости при прохождении порогового значения Яа — 0.8, в ячейке рождается одновихревой стационарный режим (рис. 2, а). При увеличении над-критичности ось симметрии основного вихря наклоняется, а в углах появляются небольшие застойные зоны, что сначала неплохо согласуется с экспериментом (рис. 3, а, [4]). Однако, имитируя в ходе численного моделирования эксперимент, а именно, медленно продвигаясь вверх по параметру, одновихревое течение при е = 0 может сохранять свою устойчивость вплоть до Яа = 1,9, после чего в системе сразу рождается автоколебательный пульсационный режим. Появление моды (2,1) при Яа » 1,0 (рис. 1, б) и ее нелинейное взаимодействие с основной модой (1,1) не приводит к перестройке течения, т.е. мода (1,1) остается доминирующей. Это противоречит опытам, в которых одновихревой стационарный режим теряет устойчивость и переходит в двухвих-ревой при Яа — 1,5 (рис. 3, в). В дополнение, чем более выраженной в эксперименте является зависимость вязкости от температуры, тем сильнее центры вихрей смещаются к низу полости [4], т.е. симметрия «верх-низ» течения нарушается, чего нет на рис. 2, б и 2, в. Еще одно противоречие связано с тем, что в расчетах, задавая начальные условия для двухвихревого движения с разной закруткой, можно на равных правах реализовать оба режима (рис. 2, б; 2, в). В то время как в эксперименте двухвихревое течение всегда характеризовалось опускным движением в центре и подъемом жидкости по бокам полости. Таким образом, можно заключить, что динамика переходов от одного режима к другому, не совсем адекватно описывается стандартными уравнениями тепловой конвекции в приближении Буссинеска. Для полноты описания было проведено численное моделирование стационарных течений в ячейке Хеле-Шоу в случае идеально теплопроводных узких вертикальных граней. Оно показало, что теплопроводность узких вертикальных граней играет немаловажную роль при формировании сценария усложнения конвективных течений в ячейке Хеле-Шоу. Однако оказалось, что пороги возникновения режимов плохо согласуются с экспериментом, и в зависимости от начальных условий может быть реализовано двухвихревое течение с обоими вариантами закрутки вихрей.

5. Результаты для неоднородной вязкости. Представим результаты расчетов для е = 0,02 1/К, что примерно соответствует воде при температуре 20 °С. При использовании граничных условий, соответствующих идеально теплопроводным узким вертикальным граням, пороговое

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

WWa fttí, . . ___

Рис. 3. Течения, наблюдавшиеся в эксперименте [4]: а) и б) - одновихревой режим при Ра = 0,9 и 1,4;

в) - двухвихревой режим при Ра = 1,5

5 10 15 20 25 5 10 15 20 25 5 10 15 20 25 5 10 15 20 25 Рис. 4. Поля функций тока стационарных конвективных течений для идеально теплопроводных узких граней в случае неоднородной вязкости; а) и б) одновихревой режим при Ра = 1,2 и Ра = 1,4; в) и г) двухвихревые режимы с «не экспериментальной» и «экспериментальной» закруткой при Ра = 1,8

Увеличение надкритичности приводит к тому, что при прохождении Ra = 1,6 система проходит точку бифуркации, в которой нижний угловой вихрь начинает вытеснять основное течение. В свою очередь, основное течение поглощает верхний угловой вихрь. В ходе такой перестройки в ячейке остаются два симметричных вихря с экспериментальной закруткой (рис. 4, г). Такой переход от одного режима к другому из-за отсутствия симметрии между угловыми вихрями теперь является единственно возможным. Второй фактор, подтверждающий роль зависимости вязкости от температуры заключается в воспроизводстве эффекта смещения центров вихрей в нижнюю часть полости (рис. 3, в и 4, г). Следует отметить, что в расчетах все же можно получить два типа двухвихревых течений с обеими закрутками (рис. 4, в; 4, г). Однако для этого требуется искусственно создать движение с противоположной закруткой и соответствующей симметрией. Это искусственно созданное движение имеет конечную амплитуду и не может быть реализовано обычными средствами в ходе эксперимента.

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

Литература

1. Путин, Г.Ф. Экспериментальное исследование надкритических конвективных движений в ячейке Хеле-Шоу / Г.Ф. Путин, Е.А. Ткачева // Изв. АН СССР. МЖГ. - 1979. - № 1. - С. 3-8.

2. Любимов, Д.В. Конвекция в ячейке Хеле-Шоу при подогреве снизу / Д.В. Любимов, Г.Ф. Путин, В.И. Чернатынский // Гидродинамика: сб. науч. тр. - Пермь, Изд-во Пермск. унта. 1977. - Вып. 10. - С. 3-14.

3. Бабушкин, И.А. Экспериментальное и теоретическое исследование переходных конвективных режимов в ячейке Хеле-Шоу / И.А. Бабушкин, В.А. Демин // Изв. РАН, МЖГ. - 2006. -№ 3. - С. 3-10.

4. Babushkin, I.A. Experimental and theoretical investigation of transitional convective flows in Hele-Shaw cell / I.A. Babushkin, V.A. Demin, D.V. Anferov // Proceeding of International Conference. «Advanced Problems in Thermal Convection». - Perm, Russia, 2004. - P. 173-178.

5. Об изменчивости одного типичного течения в ячейке Хеле-Шоу / И.А. Бабушкин, В.А. Демин, И.В. Глазкин // Изв. РАН. МЖГ. - 2009. - № 5. - С. 3-14.

6. Гершуни, Г.З. Устойчивость конвективных течений / Г.З. Гершуни, Е.М. Жуховицкий, А.А. Непомнящий. - М.: Наука, 1989. - 320 с.

7. Самарский, А.А. Теория разностных схем / А.А. Самарский. - М.: Наука, 1977. - 656 с.

8. Тарунин, Е.Л. Вычислительный эксперимент в задачах свободной конвекции / Е.Л. Тарунин. - Иркутск: изд-во Иркут. ун-та, 1990. - 228 с.

Поступила в редакцию 10 февраля 2017 г.

Bulletin of the South Ural State University Series "Mathematics. Mechanics. Physics" _2017, vol. 9, no. 2, pp. 47-54

DOI: 10.14529/mmph170206

THE EFFECT OF TEMPERATURE DEPENDENCE OF THE VISCOSITY ON STATIONARY CONVECTIVE FLOWS IN HELE-SHAW CELL

V.A. Demin, M.I. Petukhov

Perm State National Research University, Perm, Russian Federation E-mail: demin@psu.ru

The results of direct numerical simulation of stationary convective flows in a vertical Hele-Shaw cell under the uniform heating from below are presented in this paper. The calculations have been fulfilled for realistic values of the heat-transfer coefficient on vertical wide boundaries and model thermal conditions on narrow vertical walls. The approximation of plane trajectories has been applied to calculate the flows in the Hele-Shaw cell. The linear stability analysis is executed for the situation when the viscosity depends on the temperature. An analytical formula for critical values of Rayleigh number has been deduced which determine the threshold of convection in dependence on parameters of the problem. It has been shown that the numerical simulation imitating the full-scale experiment gives adequate description of the transition from one-vortex stationary flow to the two-vortex steady regime when the dependence of viscosity on the temperature is taken into account in mathematical model. The equations system of thermal convection in Boussinesq approximation was solved by the method of finite

differences at the "PGU-Tesla" supercomputer of the Research Academic Center "Parallel and Distributed Calculations" in Perm State University. The fields of stream function in vertical section have been calculated which confirm the effect of the vortices centers displacement to the bottom of the cavity.

Keywords: Hele-Shaw cell; heating from below; stationary regimes; non-homogeneity of viscosity.

References

1. Putin G.F., Tkacheva E.A. Izvestia RAN, Mekhanika Zhidkosti i Gaza, 1979, no. 1, pp. 3-8. (in Russ).

2. Lyubimov D.V., Putin G.F., Chernatynskiy V.I. Gidrodinamika, 1977, Issue 10, pp. 3-14. (in Russ).

3. Babushkin I.A., Demin V.A. Izvestia RAN, Mekhanika Zhidkosti i Gaza, 2006, no. 3, pp. 3-10. (in Russ.).

4. Babushkin I.A., Demin V.A., Anferov D.V. Experimental and theoretical investigation of transitional convective flows in Hele-Shaw cell. Proceeding of International Conference "Advanced Problems in Thermal Convection", Perm, Russia, 2004, pp. 173-178.

5. Babushkin I.A., Demin V.A., Glazkin I.V., Platonova A.N., Putin G.F. Izvestia RAN, Mekhanika Zhidkosti i Gaza, 2009, no. 5, pp. 3-14. (in Russ.).

6. Gershuni G.Z., Zhukhovitskiy E.M., Nepomnyashchiy A.A. Ustoychivost' konvektivnykh techeniy (The stability of convective flows). Moscow, Nauka Publ., 1989, 320 p. (in Russ.).

7. Samarskiy A.A Teoriya raznostnykh skhem (Theory of difference schemes). Moscow, Nauka Publ., 1977, 656 p. (in Russ.).

8. Tarunin E.L. Vychislitel'nyy eksperiment v zadachakh svobodnoy konvektsii (Computational experiment in free convection problems). Irkutsk, izd-vo Irkut. Un-ta Publ., 1990, 228 p. (in Russ.).

Received February 10, 2017