CC BY
36
Вычислительные технологии
Том 4, № 5, 1999
КОНВЕКЦИЯ И ПЕРЕНОС ТЕПЛА В ЖИДКОСТИ ПРИ ПОНИЖЕННОЙ ГРАВИТАЦИИ И УЧЕТЕ ТЕРМОКАПИЛЛЯРНЫХ ЭФФЕКТОВ *
Г. В. Глдияк, Е.А. Чеблакова Институт вычислительных технологий СО РАН Новосибирск, Россия e-mail: [email protected]
Two-dimensional flow of a viscous incompressible fluid is considered in a square laterally heated cavity with a free boundary. The effect of the gravitational and thermocapillary convection on the temperature and velocity fields distribution is investigated over a wide range of dimensionless parameter variation for a system of equations in the Boussinesque approximation. The limiting cases of dimensionless parameter variations are numerically analyzed.
Введение
Известно, что в неравномерно нагретой жидкости возникает движение: при отсутствии свободных границ — за счет тепловой (гравитационной) конвекции, при наличии свободных границ в неравномерно нагретой жидкости — за счет термокапиллярной конвекции [1, 2]. Оба процесса необходимо принимать во внимание в условиях пониженной гравитации. Представляется интересным проанализировать действие этих факторов на формирование тепловых потоков и полей скоростей, начиная с условий отсутствия гравитации и затем при увеличении ее влияния.
Будем рассматривать простейшую плоскую стационарную конвективную систему с боковым подогревом. Данная постановка задачи реализуется на практике, например, в установках по росту кристаллов [2, 3], различного типа энергетических установках. Поэтому знание структуры течения и сформировавшихся тепловых потоков в широком диапазоне изменения безразмерных параметров (числа Релея Ra, числа Прандтля Pr, числа Маран-гони Ma) представляет большой научный и практический интерес.
1. Модель
Для плоской геометрии исходные стационарные уравнения конвекции в приближении Бус-синеска в условиях невесомости для полей вихря ш, функции тока ф и температуры T в
* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований, грант №96-01-00137.
© Г. В. Гадияк, Е. А. Чеблакова, 1999.
безразмерных величинах имеют вид [4]
д , , Pr
7— (пи) + т-(vu) = f—- V2 и, (1)
дх ' дук ' Ma v ' v у
vV = и, (2) д д 1
5 (uT) + 5 (vT) = — V2 T, (3) дх ду Ma
дф дф
где п = ——, v = — — горизонтальная и вертикальная составляющие скорости, ду дх
дп 5v
и = —--—--вихрь. Данная задача характеризуется следующими параметрами: числом
ду дх
ATa L
Марангони Ma =-и числом Прандтля Pr = —. Здесь AT = Th — Tq , Th — темпера pa
ратура горячей, Tq — холодной стенки, от — температурный коэффициент поверхностного натяжения, L — характерная длина (сторона квадратной полости), ^ — коэффициент динамической вязкости, p — плотность, a — температуропроводность.
Вычисления проводились для различных чисел Прандтля и Марангони. Рассматривались также два предельных случая: Ma ^ 0 и Pr ^ ж. Если число Марангони стремится к нулю (что соответствует равенству нулю силы поверхностного натяжения), то система (1)-(3) преобразуется к виду
V2w = 0, (4)
V^ = и, (5)
V2 T = 0. (6)
Если число Прандтля стремится к бесконечности (что соответствует случаю сильно вязкой жидкости), то вместо системы (1)-(3) будем рассматривать уравнения
V2w = 0, (7)
V^ = и, (8)
д д 1
д (uT) + д (vT) = — V2 T. (9)
дх ду Ma
В случае, когда сила тяжести g отлична от нуля, исходные стационарные уравнения тепловой конвекции в приближении Буссинеска в однородном поле тяжести в переменных (ф — и) ив безразмерных величинах имеют вид [4]
д , , д , , Pr 2 Ra ■ Pr 5T
— (пи) + —(vu) = — V2 и--^тт", (10)
дх ду Ma Ma2 дх
V^ = и, (11)
i (uT ) + ly(vT ) = Ma V2 T' (12)
З _ врдЬ3(Тн - Тс) 1 др й фф
Здесь Ка =--число Релея, р = ——--температурный коэффициент
ац р дТ
расширения жидкости.
В данной работе рассматриваются следующие граничные условия (рис. 1). На одной из вертикальных стенок поддерживается постоянная температура Т = Тн = 0.5 ("горячая" стенка), на другой Т = Тс = -0.5 ("холодная" стенка). Нижняя горизонтальная стенка и свободная поверхность теплоизолированы. На твердых стенках для горизонтальной и вертикальной составляющих скорости и и V ставятся условия прилипания. На свободной
ди
поверхности компонента ц— тензора вязких напряжений должна быть равна действую-
дУ
дТ
щей на поверхность тангенциальной силе — ау——, а вертикальная составляющая скорости
дх
V — нулю. Поэтому граничные условия для рассматриваемых уравнений записываются следующим образом:
дф
ф = — = 0, Т = Тн = 0.5 при х = 0,
дх
дф
ф = — =0, Т = Тс = —0.5 при х = 1,
дх
дф дТ
ф = = 0, — = 0 при у = 0, дУ дУ
0 д2ф дТ дТ 0 .
ф = 0 ^Т = —, "5" = 0 пРи У = 1.
ду2 дх ду
Рассматривался также случай равенства нулю силы поверхностного натяжения. Это
/л\ гп\ ди
соответствует системе (4) - (6) с граничным условием ——
дУ
= 0 :
у=1
д 2ф
ф = 0, тгг = 0 при у = 1. ду2
Рис. 1. Эскиз к постановке задачи.
2. Метод решения
Постановка задачи (1) — (3), как и многих других задач вязкой несжимаемой жидкости в переменных (ф,ш) обладает следующей особенностью. Граничные условия на твердой стенке задаются только для функции тока, а не для вихря, который определен лишь внутри области согласно уравнению (2). Для преодоления этой трудности используют различные подходы, в частности, применяют приближенные граничные условия для вихря. В данной работе для вихря на твердой стенке ставилось условие Тома, которое получается из условия прилипания [5]:
2(фк+1 - фк) П/ДМ
=-Ак- + 0(Ак)- (13)
Здесь Ак — длина отрезка от граничной точки к до ближайшей к стенке точки к + 1. Использовалось также условие Вудса, которое имеет второй порядок точности. Результаты расчетов для этих двух условий близки. Однако использование приближенных граничных условий для вихря для уравнений (1)-(3) и (10)-(12) при больших числах Марангони и Релея и на подробных сетках, как следует из многих работ (например, [6-8]), связано со значительным замедлением сходимости. По этой причине в работе также применялся метод расчета граничных условий [6], который позволил улучшить сходимость в 2-4 раза, а в некоторых случаях на порядок, по сравнению с вариантом применения формулы Тома. Важным преимуществом данного метода является вычислительная устойчивость, позволяющая существенно уменьшить число итераций [6].
Идея метода состоит в том, чтобы граничное условие для вихря ставить не на границе, а внутри основной области, где вихрь определен согласно (2). Уравнение для вихря (1) решается во вспомогательной области, твердая граница которой отстоит от твердой границы основной области на один шаг сетки. Свободная граница у этих областей общая. На ней вихрь определяется следующим образом:
. ди
у= = ду
у=1
дТ дх
у=1
Уравнения для температуры и функции тока решаются в основной области. Для функции тока ставятся два краевых условия. При этом условие ф|х=о,х=1,у=о,у=1 = 0 используется при решении уравнения для ф. Полученное поле функции тока еще не удовлетворяет
дф\ п /
= 0, поэтому значения ф на границе вспомога-
ж=0, х=1, у=0
тельной области уточняются с помощью разностного аналога этого условия [6]. Используя
дф
градиентному условию ——
дп
трехточечную аппроксимацию второго порядка точности для производной —— , полу-
дп) о
чим 1
при х = 0 Ф23 = 4 фэ,, 3 = 2,...,т - 1; при х =1 фп-1, = 1 фп-2,, 3 = 2,... ,т - 1;
при у = 0 фг2 = 1 ффгэ, г = 2,...,П - 1.
При аппроксимации уравнений (1), (3) и (10), (12) для температуры и вихря применялся метод экспоненциальной подгонки [9, 10], который позволил получить более высокую
точность по сравнению с обычной аппроксимацией. В результате получилась система алгебраических уравнений с 5-точечным шаблоном. Уравнения для ш и Т не удовлетворяют условию диагонального преобладания. Известно, что при отсутствии этого условия многие эффективные методы прекращают сходиться или сходятся очень медленно. В данной работе используется модификация метода Булеева [11] и схема продольно-поперечной прогонки [12], которые позволяют находить решение и при отсутствии диагонального преобладания.
С целью улучшения сходимости для вихря ш применялась процедура релаксации [13, 14]. Она определяется следующим рекуррентным соотношением:
шп+1 = 9шп+1 + (1 — (14)
где 9 — параметр релаксации, шП — значение ш с п-й итерации после применения процедуры релаксации, шп+1 — значение ш с (п + 1)-й итерации. При численных расчетах для сетки 61 х 61 параметр 9 был порядка 0.002.
3. Анализ результатов
Вычисления проводились на прямоугольных неравномерно сгущающихся к границам области сетках размером 21 х 21, 41 х 41 и 61 х 61 в диапазоне чисел Марангони Ма от 10-3 до 104, чисел Прандтля Рг от 1 до 100 и чисел Релея И,а от 0 до 106. Исследовались также два предельных случая: Ма ^ 0 и Рг ^ ж. Результаты расчетов представлены в виде графиков и количественных характеристик течения на сетке 61 х 61. Результаты расчетов для Ма = 104 не очень достоверны из-за плохой сходимости вычислительной схемы для этого значения безразмерного параметра.
На рис. 2 показано влияние силы поверхностного натяжения на поле температур при отсутствии гравитации (И,а = 0). Если число Марангони Ма = 0, т. е. сила поверхностного натяжения равна 0 (рис. 2, а), то конвекции нет и жидкость покоится. При увеличении числа Марангони усиливается конвективное перемешивание жидкости, а при Ма = 104 течение становится вихревым. Из рис. 2, б, в, г видно, что чем больше число Марангони, тем сильнее прогибаются изолинии и тем сильнее они прижимаются к горячей и холодной стенкам. Возникают температурные пограничные слои вблизи боковых стенок, причем вблизи горячей стенки пограничный слой шире, чем вблизи холодной. Характерно также сгущение изолиний к правому верхнему углу, который соответствует области вблизи свободной границы и холодной стенки. В середине области формируется вихрь. Изолинии на рис. 2 соответствуют числу Прандтля Рг = 1 . На рис. 3 изображены изолинии температуры при Рг = ж, что соответствует сильно вязкой жидкости (система уравнений (7)-(9)). Соответствующие изолинии (см. рис. 2, 3) отличаются поведением у свободной границы у =1. В случае сильно вязкой жидкости изолинии подходят перпендикулярно к свободной поверхности и дольше остаются неизогнутыми вблизи нее в отличие от случая Рг = 1. Это объясняется тем, что силе поверхностного натяжения труднее перемешивать более вязкую жидкость. Из сравнения изолиний для Ма = 104 (рис. 2, г и 3, в) видно, что для Рг = 1 размер области вихря больше и распространен он практически на всю область. В случае Рг = ж вихрь находится ближе к свободной границе, конвективное течение во всей области слабее. Таким образом, чем больше вязкость жидкости, тем слабее там конвекция при одинаковой силе поверхностного натяжения.
На рис. 4 показано влияние гравитации на поле температур. Видно, что с увеличением числа Релея картина течения усложняется, усиливается конвективное перемешивание
Рис. 3. Изолинии температуры при Рг = ж, Иа = 0, Ма = 102 (а), 103 (б), 104 (в).
Рис. 4. Изолинии температуры при Рг =1, Ма = 102, Яа = 103 (а), 104 (б), 105 (в), 106 (г).
жидкости, появляются вихри, вблизи холодной и горячей стенок формируются узкие пограничные слои с резким перепадом температур. Сравнение соответствующих изолиний рис. 2 и 4 показывает, что сила поверхностного натяжения стремится образовать вихри в центре области и сгустить изолинии к верхней части холодной стенки, тогда как сила тяжести стремится сформировать вихри ближе к границам, оставляя середину области без вихрей, и сгустить изолинии к вертикальным стенкам. Причем при И,а = 106 образуется два узких пограничных слоя вблизи боковых стенок шириной примерно 0.07 каждый, с горизонтальным градиентом температуры. В остальной части области шириной около 0.85 изолинии располагаются более разреженно, и градиент температуры направлен вертикально вверх. Такое поведение изолиний температуры отражается на поле функции тока. Из рис. 5, г для изолиний ф (И,а = 106) видно, что вблизи боковых стенок изолинии сгущаются и образуются два узких пограничных слоя с резким изменением значения функции тока. В остальной части области изолинии располагаются более разреженно. При И,а = 105 (рис. 5, в) начинается образование вторичных вихрей в центре области.
Из табл. 1 видно, как растет максимум модуля скорости Vи2 + V2 и как меняется его положение с увеличением числа Релея. При И,а < 105 этот максимум расположен на свободной поверхности, причем с увеличением И,а он смещается к горячей стенке. При И,а = 106 максимум л/и2 + V2 расположен вблизи центра горячей стенки.
Таблица 1
Рис. 5. Изолинии функции тока при Рг = 1, Ма = 102, Иа = 103 (а), 104 (б), 105 (в), 106 (г).
Зависимость значения максимума модуля скорости и его положения от числа Релея
И,а тах V х у
102 11.871 0.6776 1
103 13.639 0.6776 1
104 26.626 0.5 1
105 73.827 0.3225 1
106 192.* 8 0.0541 0.5969
Была построена зависимость положения максимума Vи2 + V2 от числа Марангони (рис. 6). Для всех рассматриваемых чисел Ма максимум модуля скорости находится на свободной границе. С ростом Ма он смещается к правому верхнему углу (к границе холодной стенки и свободной поверхности). Для более вязкой жидкости максимум модуля скорости расположен ближе к горячей стенке, чем для менее вязкой при тех же числах Марангони. Это понятно: сила поверхностного натяжения стремится сместить вихрь к верхней границе холодной стенки, а более вязкую жидкость труднее перемешать. В промежутке изменения Ма от 103 до 104 видна осцилляция положения максимума модуля скорости (см. рис. 6), которая связана с более сложной структурой течения для больших
чисел Марангони. Для расчета течения с числами Марангони Ма > 104 требуется большая точность. С помощью предложенной схемы такой точности достичь не удалось.
Из рис. 7 для изолиний ф видно, что с увеличением числа Марангони изолинии начинают стягиваться к верхним углам, а в нижних углах области формируются вторичные вихри. С ростом Ма увеличивается размер области, которую занимает внутренний вихрь, причем модуль значения функции тока, соответствующий этому вихрю, уменьшается. У свободной границы изолинии сгущаются.
На рис. 8 представлены изолинии функции тока для сильно вязкой жидкости при различных числах Марангони. Из сравнения рис. 7 и 8 видно, что в более вязкой жидкости конвективное течение во всей области слабее, вихрь находится ближе к свободной поверхности и область его распространения меньше, в нижних углах области нет вторичных вихрей.
Сравним рис. 7, а (Рг =1, Ма = 102, Иа = 0) и 5, а (Рг = 1, Ма = 102, Иа = 103). Видно, что с появлением тепловой гравитационной конвекции меняется форма изолиний функции тока. Нижние части изолиний располагаются ближе к дну области, изолинии начинают стягиваться к двум нижним углам. Это результат влияния силы тяжести. Чем
Рис. 6. Зависимость положения максимума модуля скорости от числа Марангони.
Рис. 7. Изолинии функции тока при Рг =1, Яа = 0, Ма = 102 (а), 103 (б), 104 (в).
больше число Релея, тем сильнее влияние силы тяжести по сравнению с влиянием силы поверхностного натяжения на характер течения. Так, при И,а = 103, Ма = 102 (см. рис. 5, а) изолинии также стягиваются к правому верхнему углу, т. е. термокапиллярная конвекция
оказывает сильное воздействие на форму изолиний. С дальнейшим увеличением Иа гравитационная конвекция становится преобладающей. Начиная с Иа =105 происходит образование вторичных вихрей, причем один стремится в левый верхний угол области, к более нагретой стенке, другой — в правый нижний, к холодной стенке (см. рис. 5, б - г).
Были выполнены расчеты задачи конвекции в условиях невесомости при нулевой силе
поверхностного натяжения (уравнения (4) - (6) с нулевым граничным условием —— =
ду у=1
0). В этом случае на жидкость не действуют никакие силы, конвекция не возникает и все профили функции тока, скоростей и вихря вдоль различных сечений области практически нулевые. Изолинии температуры в этом случае являются прямыми линиями (см. рис. 2, а). Это говорит о том, что тепло передается только за счет теплопроводности.
Были проведены также расчеты со следующими параметрами: Ма = 1; 0.1; 0.01; 0.001; Рг =1. С уменьшением числа Марангони профили функции тока и скоростей монотонно стремятся к нулевым значениям, к профилям предельного случая Ма = 0. Изолинии температуры с уменьшением числа Ма также монотонно приближаются к изолиниям случая Ма = 0.
Наряду с этим рассматривался другой предельный случай: Рг ^ го для различных чисел Марангони в условиях отсутствия силы тяжести (уравнения (7)-(9)). Результаты расчетов показали, что для любого числа Марангони в промежутке Ма от 0 до 10 течение не зависит от числа Прандтля (при Рг от 1 до го). Для Ма = 102 и Рг от 10 до го соответствующие профили течения совпадают (см. рис. 3, а (изолинии температуры), рис. 8, а (изолинии ф)), но отличаются от профилей с числом Прандтля Рг = 1 (см. рис. 2, б, 7, а). То же самое происходит при Ма = 103 (см. рис. 2, в, 7, б; 3, б; 8, б). Для Ма = 104 и Рг = 100, го профили близки между собой и отличаются от профилей течения при Рг = 10.
Для контроля теплового баланса вычислялись интегральные числа Нуссельта вдоль горячей и холодной стенок и вдоль сечения х = 1/2. Локальное число Нуссельта в точке с координатами (х, у) и интегральное число Нуссельта МиХ0, соответствующее сечению х = х0, определяются равенствами
дТ
Ш(х,у) = Ма • иТ - —,
дх
N"^0 = / Nu(xo,y)dy.
Зо
В табл. 2 приведены результаты вычислений задачи конвекции в области со свободной границей и боковым подогревом в условиях невесомости для Рг =1 и различных чисел Марангони. Здесь и(1/2,1) — значение скорости и в центре свободной поверхности, N"0, N"1/2, N"1 — соответственно интегральные числа Нуссельта вдоль сечений х = 0, х = 1/2, х = 1. Полученные результаты сравниваются с результатами [15]. В работе [15] для решения уравнений (1) - (3) применялась схема повышенного порядка точности. Из таблицы видно, что используемый в данной работе метод позволяет получить с приемлемой точностью интегральные характеристики теплоотдачи (N"0, N"1/2, N"1). Однако локальные характеристики (N"(0,1), N"(1,1)) определяются с существенными погрешностями, особенно при Ма = 104.
Таблица 2
Результаты вычислений для Pr = 1
Ма Сетка и(1/2,1) N"0 N"1/2 N"(0,1) N"(1,1)
102 21 х 21 1.0326 -1) 1.1322 1.0999 1.0396 0.7972 1.7816
41 х 41 1.0711 -1) 1.1010 1.0977 1.0867 0.7402 1.9874
61 х 61 1.0671 -1) 1.0945 1.0953 1.0924 0.7313 1.9905
[15] 61 х 61 1.0869 -1) 1.0962 1.0962 1.0962 0.7301 2.028
103 21 х 21 5.2221 -2) 2.1800 1.9976 1.6900 1.3444 7.3224
41 х 41 5.2032 -2) 1.9858 1.9742 1.8983 1.0268 10.6884
61 х 61 5.1404 -2) 1.9437 1.9528 1.9297 0.9707 11.438
[15] 61 х 61 5.0018 -2) 1.9258 1.9258 1.9258 0.9550 11.75
104 21 х 21 2.7263 -2) 5.4250 3.9362 2.6086 5.1195 12.6693
41 х 41 3.1541 -2) 4.5622 4.3744 4.1286 3.1227 38.5290
61 х 61 3.1562 -2) 4.4459 4.4467 4.3524 2.5091 60.471
[15] 61 х 61 3.0381 -2) 4.3621 4.3621 4.3654 2.2334 77.09
Были также построены зависимости интегральных чисел Нуссельта от чисел Пранд-тля, Марангони и Релея (соответственно рис. 9, 10, 11). Из зависимости ^ от Рг (см. рис. 9) для различных значений Ма видно, что чем больше число Марангони для одного и того же числа Прандтля (т.е. чем больше термокапиллярная конвекция), тем сильнее теплоотдача. Для Ма от 0 до 102 коэффициент теплоотдачи не меняется с ростом Рг. При Ма = 103; 104 и начиная с Рг =10 он выходит на некоторый постоянный закон, причем для каждого числа Ма этот закон свой. Кроме того, в промежутке Ма от 103 до 104 наблюдается скачок числа Нуссельта, особенно для случая Рг =1. Причиной этого, по-видимому, является возникновение сложной структуры течения и потеря им устойчивости [15], вследствие чего резко увеличивается теплоотдача. Такой скачок теплоотдачи наблюдается и на рис. 10, характеризующем зависимость Nu от Ма для различных Рг. Число Нуссельта увеличивается с ростом числа Релея, характеризующего интенсивность гравитационной конвекции (см. рис. 11). В промежутке И,а от 104 до 106 также происходит скачок Nu из-за возникновения сложной структуры течения. На рис. 10, 11 треугольниками показаны функции, аппроксимирующие зависимость числа Нуссельта вдоль горячей стенки от числа Релея в промежутке И,а от 105 до 106 и зависимость числа Нуссельта вдоль горячей стенки от числа Марангони в промежутке Ма от 103 до 104. Эти функции получе-
Рис. 9. Зависимость интегрального числа Нуссельта вдоль горячей стенки от числа Прандтля.
Рис. 10. Зависимость интегрального числа Нуссельта вдоль горячей стенки от числа Марангони (g(Ma) = = 1.481 + 7.15 • 10-4 • Ma -- 3.277 • 10-8 • Ma2).
Рис. 11. Зависимость интегральных чисел Нуссельта вдоль горячей и холодной стенок от числа Релея (0(Бя) = 4.122+9.507• 10-6 х х Ra - 4.296 • 10-12 • Ra2).
ны методом наименьших квадратов и являются линейными комбинациями функций вида
#(Ма) = 1.481 + 7.15 • 10-4 • Ma - 3.277 • 10-8 • Ma2 -в случае зависимости числа Нуссельта от числа Марангони и
0(Иа) = 4.122 + 9.507 • 10-6 • Иа - 4.296 • 10-12 • Иа2 -в случае зависимости числа Нуссельта от числа Релея.
4. Заключение
Исследовано влияние гравитационной и термокапиллярной конвекции на поле температур и скоростей. Расчеты показали, что начиная с Иа = 104 гравитационная конвекция становится преобладающей при Ма < 102.
С ростом гравитационных сил картина течения усложняется, усиливается конвективное перемешивание жидкости, вблизи холодной и горячей стенок образуются узкие пограничные слои с резким перепадом температур и скоростей. Начиная с Иа = 105 происходит образование вторичных вихрей.
Конвективное течение во всей области слабее для более вязкой жидкости, вихри находятся ближе к свободной поверхности, и область их распространения меньше.
Рассматривались два предельных случая: Ма ^ 0 и Рг ^ то в условиях отсутствия силы тяжести. Расчеты показали, что с уменьшением числа Марангони профили течения монотонно приближаются к профилям течения для случая Ма = 0. Также были определены диапазоны изменения безразмерных параметров, в которых течение не меняется.
С помощью метода наименьших квадратов были построены функции, аппроксимирующие зависимости числа Нуссельта вдоль горячей стенки от числа Релея в промежутке Иа от 105 до 106 и от числа Марангони в промежутке Ма от 103 до 104.
Авторы выражают благодарность В. С. Бердникову, В. А. Гапонову, Г. Б. Чеблакову, А. В. Двуреченскому за полезные советы и обсуждения в ходе работы.
Список литературы
[1] БЕРДНИКОВ В. С., Кирдяшкин А. Г. Структура свободноконвективных течений в горизонтальном слое жидкости при различных граничных условиях. В "Структура пристенного пограничного слоя". Новосибирск, 1978, 5-45.
[2] ПОЛЕЖАЕВ В. И. Гидродинамика, тепло- и массообмен при росте кристаллов. В Ав-дуевский В. С., Полежаев В. И. "Гидромеханика и тепломассообмен при получении материалов". Наука, М., 1990.
[3] Дубовик К. Г., Никитин С. А., ПОЛЕЖАЕВ В. И. Гидродинамические эффекты температурного и концентрационного расслоения. В "Проблемы вязких течений. Изд. ИТПМ СО РАН, Новосибирск, 1981, 55-64.
[4] ЛАНДАУ Л. Д., Лифшиц Е. М. Теоретическая физика. Том 6. Гидродинамика. Наука, М., 1986.
[5] РОУЧ П. Вычислительная гидродинамика. Мир, М., 1980.
[6] ПОЛЕЖАЕВ В. И., ГРЯЗНОВ В. Л. Метод расчета граничных условий для уравнений Навье — Стокса в переменных "вихрь, функция тока". Докл. АН СССР, 219, №2, 1974, 301-304.
[7] КУСКОВА Т. В. Численное исследование двумерных течений вязкой несжимаемой жидкости. В "Некоторые применения метода сеток в газовой динамике". Вып. 3, 1971.
[8] ПОЛЕЖАЕВ В. И., Вальциферов Ю. В. Численное исследование нестационарной тепловой конвекции в цилиндрическом сосуде при боковом подводе тепла. Там же.
[9] CHEBLAKOVA E.A., Gadiyak G.V. Convective flow in an upright enclosed square cavity: a comparison exercise. IJCFD, 10, No. 2, 1998, 139-149.
[10] МАРЧУК Г. И. Методы вычислительной математики. Наука, Новосибирск, 1973.
[11] SCHNEIDER G. E., Zedan M. A modified strongly implicit procedure for the numerical solution of field problems. Numerical Heat Transfer, 4, 1981, 1-19.
[12] ЯНЕНКО Н. Н. Метод дробных шагов решения многомерных задач математической физики. Наука, Новосибирск, 1967.
[13] Тарунин Е. Л. Анализ аппроксимационных формул для вихря скорости на твердой границе. Уч. зап. Пермск. гос. пед. ин-та. Гидродинамика. №152, вып. 9,1976,167-179.
[14] МызниКОВА Б. И., Тарунин Е.Л. О граничных условиях для вихря скорости в задачах динамики вязкой жидкости. В "Конвективные течения и гидродинамическая устойчивость". Изд. УНЦ АН СССР, Свердловск, 1979, 90-101.
[15] Бердников В. С., ГАПОНОВ В. А. Тепловая гравитационно-капиллярная конвекция в прямоугольных полостях. Процессы тепломассопереноса и рост монокристаллов и тонкопленочных структур. Тр. II Российского симпозиума (HT&CG'97). Обнинск, 1998.
[16] Самарский А. А., Николаев е. С. Методы решения сеточных уравнений. Наука, М., 1978.
[17] Годунов С. К., Рябенький В. С. Разностные схемы. Наука, М., 1973.
[18] Пейре Р., Тейлор Т.Д. Вычислительные методы в задачах механики жидкости. Гидрометеоиздат, Л., 1986.
[19] Полежаев В. И., Бунэ А. В., Верезуб Н. А. и др. Математическое моделирование конвективного тепломассообмена на основе уравнений Навье — Стокса. Наука, М., 1987.
Поступила в редакцию 8 октября 1998 г., в переработанном виде 27 января 1999 г.
