Моделирование конвекции в областях со свободными границами Текст научной статьи по специальности «Физика»

Вычислительные технологии

Том 5, № 6, 2000

МОДЕЛИРОВАНИЕ КОНВЕКЦИИ В ОБЛАСТЯХ СО СВОБОДНЫМИ ГРАНИЦАМИ *

Е. А. ЧЕБЛАКОВА Институт вычислительных технологий СО РАН Новосибирск, Россия e-mail: lena@ict.nsc.ru

The numerical algorithm for the plain heat convection problem for viscous incompressible liquid with free boundaries is presented. It is based on the transformation of domains with the curvilinear moving boundaries to the rectangular domains. To test the algorithm the following problems were solved numerically: convection problem for the liquid with a free boundary liquid/gas and the Stephan problem for the solid phase fusion.

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

В работе [6] численно решена стационарная задача Стефана в области со свободной границей жидкость — газ в плоском случае, где осуществляется преобразование одной пространственной координаты для выпрямления свободной границы жидкость — газ, а задача Стефана решается методом сквозного счета. Представляет интерес построить модель, в которой обе пространственные координаты преобразуются для выпрямления как границы фазового перехода, так и свободной границы жидкость — газ. В настоящей работе выполнено такое преобразование координат и получены уравнения в новых переменных. На их основе построен вычислительный алгоритм решения задач тепловой конвекции вязкой несжимаемой жидкости со свободными границами и численно решены следующие тестовые задачи, являющиеся составными частями задачи моделирования роста кристаллов методом бестигельной зонной плавки, а именно:

— задача расчета конвективных течений жидкости со свободной границей жидкость — газ;

— задача Стефана о фазовом переходе в области с фиксированными границами жидкость — газ.

Конвективные течения в этих задачах вызваны эффектами плавучести (конвекция Грасгофа) и термокапиллярности (конвекция Марангони). В первой из них исследуется влияние термокапиллярной конвекции на картину течения при боковом подогреве в

* Работа выполнена при финансовой поддержке федеральной целевой программы "Интеграция", проект №274.

© Е. А. Чеблакова, 2000.

условиях слабой гравитации, во второй — процесс плавления твердого вещества в двухслойной среде при подогреве снизу. Задача Стефана решается методом явного выделения границ раздела фаз.

Эти задачи нелинейные: неизвестная граница определяется по неизвестному решению и, наоборот, решение ищется в неизвестной заранее области.

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

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

А.

В.

С.

Са = x//, Па = у/дх.

Св = я//, (1)

у - дх . , /0\

Пв =--+ 1. (2)

д2 - дх

Сс = я//, (3)

у - д2 . 0

Пс = 7-- + 2. (4)

7 — д 2

х о 1 У о 1С

Рис. 1. Метод бестигельной зонной плавки. Вертикальный вариант (преобразование координат) (а) и горизонтальный (б).

Замечание 1. Предполагаем, что кривые д\, д2, f таковы, что с помощью указанных выше преобразований их можно сделать координатными линиями.

Запишем исходные уравнения тепловой конвекции в однородном поле тяжести в приближении Обербека — Буссинеска для жидкой фазы (область В) в переменных функция тока — вихрь и безразмерном виде [7]

ди д , ч д , ч Аи Ог дт

ж + дХ (ии) + ду м = ке - и? дХ' (5)

дт д д Ат

эТ + 55(ит) + зЦ(ит) = ил' (6)

А^ = и, (7)

д-0 д-0

где и = -г—, V = — — горизонтальная и вертикальная составляющие скорости; ду дх

ди дv

и = —--— — вихрь скорости; Ие, Ог, Рг — числа Рейнольдса, Грасгофа и Прандтля

ду дх соответственно.

В областях твердой фазы выполняется уравнение теплопроводности

дт Ат

(8)

дТ ИеРг

После преобразования координат уравнения примут вид: Область В

ди 1 / д / ди ^ ди ^ д-0 \ - ' ' Вп — + В12— — Иеи^ +

дТ Ref 72 у д£ \ д£ дп дп

д ди л ди .

+ дП (В12 дё + В22 дП + Кеид^' 1 +

+ ди д^, + (п — 1)(д2г— д1 г) — (дт + ВвдТ) (9)

дп д2 — д1 Ие2/ \ дё В дп

дт 1 /д ( дт дт д^Л

1 1 Вц— + В12дт — ИеРгт-^ +

дТ ИеРг/ 72 \д£ \ д£ дп дп

+ - (в12дт + В22дт + ИеРгт^^ + дтд;г + (п — 1)(д2г — д1 г), (10)

дп \ дё дп дё / / дп д2 — д1 '

1 1(вц Ш + ВИ ^ + £(вв + В22 = и. (11)

11 12 12 22 /Л \дё \ дё дп/ дп \ дё дп

Здесь и далее индексы у переменных п опущены; 72 = д2 — д1 + п(д2п — д1 п) — д2п + 2д1 п;

п(д2? — д;1 ?) — д2? + 2д1 ^ _ ё 0 _ 7

ВБ = ^^-^. р = ^/п. Вц = ^(1 + р2);

В12 = — р + ВцВв; В22 =

/ 1 + В?2 72 1+ р2 .

Область С

дТ = 1 (£ (с„дТ + С,2дТ) + « (СидТ+С22дТ^ + , (12)

д£ КеРг/<02 у дС \ д£ дп / дп \ д£ дп // дп Ь — д2

где

, д2$ (з — п) С ,

<о,2 = Ь — д2 + д2 п (3 — п); Во =--;-; Р = -¡— ;

<0,2 ¿0,2

Си = ^ (1+ Р2); С12 = —р + Сп Во; С22 = / .

/ <¿0,2 1 + Р2

Уравнения тепловой конвекции и уравнение теплопроводности в новых переменных £, П записаны в консервативном виде.

Замечание 2. Поскольку функции д1 и д2 являются функциями не только пространственных координат, но и времени, переменная п также зависит от времени. В связи с этим для каждой из областей А, В и С в правой части уравнений для температуры и

дТ дп ди дп

вихря появляются слагаемые ——— и ——— соответственно.

дп д^ дп д^

Существует еще так называемый горизонтальный вариант бестигельной зонной плавки [8] (рис. 1, б), для которого автором настоящей работы сделано преобразование координат и получены уравнения тепловой конвекции в новых переменных, аналогичные уравнениям (9)-(12). На основе этих уравнений построена математическая модель для решения задач тепловой конвекции с подвижными границами. Для тестирования модели решались следующие две задачи.

2. Задача 1. Исследование конвекции и переноса тепла в жидкости с подвижной свободной границей жидкость — газ

Физико-математическая модель. Пусть вязкая несжимаемая жидкость с плотностью р, кинематической вязкостью V и коэффициентом поверхностного натяжения ст(Т) занимает область В (рис. 1, б). Предположим, что коэффициент ст(Т) зависит от температуры линейно:

1 аст

ст(Т) = сто (1 — (Т — То)), сто = ст(Т0), стт = — ~т=, .

сто а! т=т0

Рассмотрим частный случай, когда д1 = 1, д2 = 2 и у = 0 — твердые непроницаемые стенки, у = /(ж) — свободная граница. Движение жидкости и теплоообмен описываются уравнениями тепловой конвекции (5)-(7).

Найдем квазистационарное решение, для чего на каждом шаге по времени считаем нормальную компоненту скорости на свободной границе vn|y=f(х) = 0. Для касательной составляющей скорости у3 на границе /(ж) выполняется соотношение [9]

ду3 ду3 2 д2у3

аГ + 1Й = И И? + г' (13)

^ ди — ^ + *( 1Л — сттА! ^ + / (&г — с).

Ие дп дв Ие дв \Я \ ст0 )) + /

,2 ' X

Здесь Р0 — внешнее давление; Са — капиллярное число; К — радиус кривизны свободной границы; О — число Галилея.

Перейдем к системе координат, в которой криволинейная граница f (ж) является координатной линией п =1 прямоугольной сетки (рис. 1, б), с помощью следующих формул, аналогичных формулам (1), (2):

lf с ж - 91 , т П = У/f, С =-+ 192 - 91

В этой системе координат уравнения тепловой конвекции имеют вид ди 1 / 5 / ди л ди д-\ д / ди л ди дф\ .

д£ Ref J2 у дС \ дС дп дп / дп \ дС дп дС дид'и + (С - 1)(д2г - 91 г) Сг /дТ п дТ ^ дТ\

+ дС 92 - 91 дС /Ч дп + ^ дС )

дТ 1 / д ( дТ дТ дф\

1 1 Б11 — + Б12 т- - ИеРгТ-^ +

д£ RePгf J2 \ дС \ дС дп дп

СбидТ + Б22£ + ИеРгТ?^|| + дТ 9!• + (С - - 9и)

дп V дС дп дС / / дС 92 - 91

1 ^4 (бц ^ + БИд-I + 4- (БИд- + Б22^Н = и,

где

^ Т \ О >= \ "^11 1 12 о I 1 О \ "^12 г\ >= 1 22 о

/Л \дС V дС дп / дп \ дС дп

т = 9 9 + С(9 9 ) 9 + 29 . С = С(92п - 91П) - 92п + 291п . ^2 = 92 - 91 + С(925 - 915) - 92? + 291С. =--г-.

Р = уЛ; #22 = ^ (1+ р2); Б12 = -р + Б22СВ; #11 = 4. J2 / J2 1 + Р2

Пусть в начальный момент времени жидкость находится в состоянии покоя. Начальные условия для температуры заданы в виде Т(ж,у, 0) = 1.5 - х при 0 < у < 1. Граничные условия на твердых стенках следующие: — = 0, д—/дп = 0,

Т = 0.5 при ж =1, Т = -0.5 при ж = 2, Т =1.5 - ж при у = 0.

Для связи функции тока и вихря на твердых стенках используется условие Тома [10], которое после преобразования координат в новых переменных С, п для границы п = 0 примет вид

'#22\ 2-4)2

и

М/ ¿,1 к

Б11 \ 2-2,7 _ /Бц\ 2-м-1,7

для границ £ = 1, £ = 2. Здесь Бц, В22 — коэффициенты квадратичной формы; / 32 — якобиан преобразования; Н£1, Нщ — шаги сетки.

На свободной границе условия для ф и ш в переменных х, у заданы соотношениями [6, 9]:

дф = 13 =2 Ма дТ

ду у/1 + /2' Я 3 Ие дэ

(у3 — решение уравнения (13); Ма — число Марангони). Считаем, что 13(1) = 13(2) = 0, Ро = 0.

Для температуры на свободной границе задано условие дТ/дп = 0. Положение свободной границы определяется использованием кинематического условия [9]

д/ дф

Ж + 1 + /2 = »■

Следуя [6], запишем граничные условия для функции тока, вихря и кинематическое условие на свободной границе в переменных £, п:

дф = (14) дп Б22 ' 1 '

2 Ма дТ

Ш = У3--/ я?! (15)

Я ие^/г/ д£

д/ , дф п

аь + д£ = 0- (16)

Уравнение для касательной составляющей скорости 13 на свободной границе (13) в переменных £, п преобразуется к виду

£+131 = ие | (уТ/ (17)

' = и1е (В121 + В22Щ) - + £| (Я (1 - ^т)) + (от - о.

Условие для температуры на свободной границе в новых переменных дТ/дп = 0 запишем так:

В12дт + Б22дт ) = 0.

I I "*—' 12 о >= 1 "*—' 22 г\

^Т/ V д£ дп

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

дф

дП У=/(X)

и уравнение (13) для касательной компоненты скорости на свободной границе.

Опишем общий алгоритм решения задачи. В прямоугольнике, соответствующем преобразованной области В, строится равномерная прямоугольная сетка. Пусть в момент времени Ьк известны значения Тк, шк, фк, /к, . Для перехода к моменту времени Ьк+1 = Ьк + т совершаем следующие действия:

1. Решаем уравнение (16) для расчета нового положения свободной границы.

2. По полученным значениям /к+1 определяем матрицы коэффициентов В11(г,]), В12(г,]), В22(%,]) в момент времени Ьк+1.

3. Решая уравнение (17) для у3, находим условия (14), (15) для ф и ш на свободной границе.

На каждой итерации последовательно решаем уравнения для температуры, вихря и функции тока по схеме стабилизирующей поправки [11], порядок аппроксимации равен 0(т + к2).

4. Определяем поле температуры, функцию тока берем с предыдущей итерации.

5. Решаем уравнение для вихря со значениями температуры, взятыми с предыдущего временного шага Ьк, и значениями ф с предыдущей итерации.

6. Вычисляем значение функции тока, причем уравнение для ф решается до установ-

ления, т.е. пока не будет достигнуто неравенство шах.

г,з

фр+1 - фр .

фР+1

< £ (р — номер ите-

рации, £ — заранее заданная точность, которая при численных расчетах равнялась Ю-5). После этого итерации 4-6 повторяются до достижения заданной точности для всех функций Т, ш, ф. На этом переход на новый временной слой Ьк+1 завершен. Далее процесс повторяется, пока не будет получено стационарное решение, которое считается найден/ _ /к з г з г

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

< £ для всех г = 1, ... , N (ф

заданное число шагов по времени). Также при достижении стационарного решения нормальная составляющая скорости уп на границе / (х) (или на координатной линии п = 1) должна быть близка к нулю.

Результаты расчетов. На рис. 2 представлены результаты расчетов для установившегося течения, проведенных на сетке 21 х 21 при И,е = 1, Ог = 0,0 = 0, Рг = 0.73, Са-1 = 1, иллюстрирующие эффекты термокапиллярности (конвекция Марангони). Видно, что в области образуется один вихрь, который "прижимается" к свободной границе. Вблизи нее изолинии функции тока сгущаются, появляется пограничный слой. С увеличением числа Марангони деформация свободной границы усиливается. В целом конвективное течение во всей области достаточно слабое в отличие от течений, вызванных эффектами плавучести [2]. Качественная картина течения и количественные характеристики (деформация свободной границы) хорошо согласуются с результатами, полученными в [6, 12]. В [12] задача со свободной границей решается в переменных и, V, р, Т методом конечных элементов. Основная особенность метода расчета, использованного в настоящей работе, в отличие от предложенного в [6], состоит в том, что мы рассматриваем более сложную геометрию области и строим общую модель.

3. Задача 2. Решение двумерной задачи Стефана о фазовом переходе

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

Рис. 2. Изолинии функции тока (a) и температуры (б) при боковом подогреве жидкости с верхней свободной границей в условиях невесомости. Ma = 1 (слева), Ma = 2 (справа).

граница раздела фаз, gi = 1, f = 1. Пусть в начальный момент времени жидкость находится в состоянии покоя, а начальное положение границы раздела фаз д2 = 2. Исследуем численно процесс плавления твердого вещества при условии, что расплав подогревается снизу.

В области B имеет место конвективное движение жидкости, которое описывается системой уравнений тепловой конвекции в приближении Обербека — Буссинеска (5)-(7). В области твердой фазы выполняется уравнение теплопроводности (8). Искомые функции — вихрь скорости, функция тока, температура и граница раздела фаз д2.

Последовательно берем начальные распределения температуры в области B:

nx

TB(x,y, 0) = 1 - y - 0.01 sin ny cos —, (18)

nx

Tb (x, y, 0) = 1 - y + 0.01 sin ny cos—, (19)

TB(x, y, 0) = 1 - y + 0.01 sin ny cos nx. (20)

В области C температура в начальный момент времени распределяется по линейному закону — от T = 0 на нижней границе y = g2(x) до T = -0.5 на границе y = 3. Граничные условия для температуры следующие: T =1 при y =1;

T = 0 при y = g2(x) (температура кристаллизации или плавления); T = -0.5 при y = 3.

Вертикальные границы x = 0 и x = 1 теплоизолированы. Для скоростей (на границе области жидкой фазы) ставятся условия скольжения: нормальная составляющая скоро-

сти на границе равна нулю. Для переменных "вихрь, функция тока" граничные условия таковы: ф = 0, и = 0.

На границе фазового перехода д2 выполнены условия сопряжения [8]:

Tb(x, g2(t, x)) = Tc(x, g2(t, x)) = Т*.

Здесь Tb и Tc — температура в областях B и C соответственно; Т* — температура плавления. На границе д2 выполняется условие Стефана [8]:

-paVn = limr<r^Н^ - Шп (Н^ + PcVnT*, (21)

где положительная постоянная а — скрытая удельная теплота плавления; Н — коэффициент теплопроводности; Vn — скорость перемещения границы фазового перехода в направлении нормали к ней; c — теплоемкость. В безразмерных величинах и с учетом того, что Т* = 0, условие Стефана (21) примет вид

V к дТв к dTc

Vn = k^ —--Kc-^t—

on on

(кв и kc — безразмерные коэффициенты теплопроводности в областях B и C соответственно) .

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

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

Метод расчета. Идея метода решения основана на декомпозиции сложной системы (многослойной среды) на более простые элементы (модули) при определенных выше граничных условиях [4]. В каждом модуле численный расчет проводится независимо от других модулей, а условия на общих границах объединяют их в исходную систему.

По формулам (1)-(4) осуществим преобразование координат так, чтобы криволинейная граница раздела фаз д2 стала координатной линией прямоугольной системы координат. Уравнения тепловой конвекции для области B и теплопроводности для твердой фазы в новых переменных ц имеют вид (9)-(11) и (12) соответственно.

Условие Стефана в переменных ц перепишем так:

Vn = TTTBsMB-f + B22f) - Ч-f + О). (22)

По формуле (22) находим нормальную составляющую скорости движения границы фазового перехода д2. Пусть — положение границы д2 в момент времени tk, д2к+1 — положение д2 в момент времени tk+1 (рис. 3), т = tk+1 — tk, тогда

дк+1 = дк + TVn/ cos а, (23)

1 11

cos а = . -= —. = = —. . (24)

Vl + tan2 а V1 + (д2 x)2 x/l+B^

Рис. 3. Схема движения границы фазового перехода $2-Кривые 1, 2 — положение $2 в момент времени и Ьк соответственно.

Подставляя (22) и (24) в (23), получим

<Г = ак + тУп= ,,к + Т (кВ (В12дГ^ + В22дГ^) - кс (С1адГ^ + С22) -(25)

По формуле (25) находим положение границы д2 в момент времени Итерационный процесс решения задачи строится следующим образом. Пусть в момент времени ¿к известны значения Тк, ск, фк, д2к.

1. По условию Стефана (22) определяем скорость УП движения границы д2 в направлении нормали к ней.

2. Находим новое положение границы раздела фаз дк+1.

3. Определяем матрицы коэффициентов Вц, В12, В22.

4. Вычисляем температурное поле Тс в области твердой фазы С.

Итерационный процесс осуществляется так же, как в задаче 1. Он заканчивается, когда достигается заданная точность для всех функций ТВ, с, ф. На этом переход на новый временной слой ¿к+1 завершен. Далее процесс повторяется начиная с п. 1, пока не будет получено стационарное решение, которое считается найденным, если выполняются условия

тах

Тк+1 _ тк

Т г,] Т г,]

Т к+1

< е,

тах

, ,к+1 / к

си,

к+1 г,]

< е,

тах

г,]

дк+1 д2 г

д2к

дк+1 д2 г

< е,

где е — заранее заданная точность, которая при численных расчетах равнялась 10-5.

Основная сложность задачи заключается в том, что коэффициенты уравнений (9) -(12) зависят от положения границы д2, которая, в свою очередь, зависит от температуры, а значит, и от всех искомых функций. Таким образом, наряду с нелинейностью по ф имеем нелинейность по д2. Для разрешения этих нелинейностей используется описанный выше итерационный процесс. За начальное приближение Уп берется значение с нижнего временного слоя.

Каждое из уравнений для Тс, ТВ, ф, сс решается методом стабилизирующей поправки. В результате получается система конечно-разностных уравнений, аппроксимирующая исходные дифференциальные уравнения с точностью 0(т + ). Эта система уравнений имеет трехдиагональную структуру и эффективно решается методом прогонки [13].

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

Рис. 4. Изолинии функции тока (а) и температуры (б) при решении задачи Стефана

с начальным распределением температуры:

пх

I — Tb (х, у, 0) = 1 — у — 0.01 sin пу cos —;

пх

II — Tb (х,у, 0) = 1 — у + 0.01 sin пу cos —;

III — Tb (x, у, 0) = 1 — у + 0.01 sin пу cos пх.

Gr = 104, Pr = Re = 1.

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

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

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

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

Автор выражает глубокую благодарность А. С. Овчаровой, Г. В. Гадияку и Г. Г. Черных за постоянное внимание к работе и ценные советы.

Список литературы

[1] САМАРСКИЙ А. А., МОИСЕЕНКО Б. Д. Экономичная схема сквозного счета для многомерной задачи Стефана // Журнал вычисл. математики и мат. физики. 1965. Т. 5, №5. С. 816-827.

[2] ТАРУНИН Е. Л. Вычислительный эксперимент в задачах свободной конвекции. Иркутск: Изд-во Иркут. ун-та, 1990.

[3] ВАБИЩЕВИЧ П. Н., ИлИЕВ О. П. Численное решение сопряженных задач тепло- и массопереноса с учетом фазового перехода // Дифференц. уравнения. 1987. Т. 23, №7. С. 1127-1132.

[4] ОВЧАРОВА А. С. Метод решения двумерной многофронтовой задачи Стефана // ПМТФ. 1995. Т. 36, №4. С. 110-119.

[5] БУДАК Б.М., УСПЕНСКИЙ А. Б. Разностный метод с выпрямлением фронтов для решения задач типа Стефана // Журнал вычисл. математики и мат. физики. 1969. Т. 9, №6. С. 1299-1315.

[6] ОВЧАРОВА А. С. Численное решение стационарной задачи Стефана в области со свободной границей // Вычисл. технологии. 1999. Т. 4, №1. С. 88-99.

[7] ГЕРШУНИ Г. З., ЖУХОВИЦКИЙ Е.М., Непомнящий А. А. Устойчивость конвективных течений. М.: Наука, 1989.

[8] МеЙРМАНОВ А. М. Задача Стефана. Новосибирск: Наука, 1986.

[9] ПУхНАЧЕВ В. В. Движение вязкой жидкости со свободными границами: Учеб. пособие. Новосибирск: Изд-во Новосиб. ун-та, 1989.

[10] РОУЧ П. Вычислительная гидродинамика. М.: Мир, 1980.

[11] ЯНЕНКО Н. Н. Метод дробных шагов решения многомерных задач математической физики. Новосибирск: Наука, 1967.

[12] Chippada S., Jue T.C., Ramaswamy D. Finite element simulation of combined buoyancy and thermocapillary driven convection in open cavities // Int. J. Num. Meth. Eng. 1995. V. 38. P. 335-351.

[13] Самарский А. А., НИКОЛАЕВ Е. С. Методы решения сеточных уравнений. М.: Наука, 1978.

Поступила в редакцию 25 июля 2000 г.