Численное решение стационарной задачи Стефана в области со свободной границей Текст научной статьи по специальности «Физика»

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

Том 4, № 1, 1999

ЧИСЛЕННОЕ РЕШЕНИЕ СТАЦИОНАРНОЙ ЗАДАЧИ СТЕФАНА В ОБЛАСТИ СО СВОБОДНОЙ ГРАНИЦЕЙ*

А. С. Овчарова Институт гидродинамики им. М.А. Лаврентьева СО РАН

Новосибирск, Россия e-mail: [email protected]

The numerical algorithm for plain problem of heat and mass transfer based on the combination of two methods is presented in this paper. In order to determine the phase boundary between a liquid and a solid, the standard approach with "smoothing"of coefficients is used. Computation of the flow in the time-dependent domain is based on the method of fictive domains. In order to find the position of the liquid free surface, a new method, which takes into account the peculiarities of the boundary conditions on the free surface, has been developed. The results of calculations of some sample problems are presented.

Исследования конвективных течений, вызванных эффектами плавучести (конвекция Грасгофа) и термокапиллярности (конвекция Марангони), получили большое развитие из-за их важности в производстве кристаллов путем бестигельной зонной плавки (БЗП). В большинстве ситуаций термокапиллярная конвекция достаточно мала в сравнении с естественной конвекцией и часто опускается, но ее вклад резко возрастает при слабой гравитации или при расчетах, проводимых для жидкого моста БЗП. Теплообмен и возникающие при этом конвективные течения в расплаве не только влияют на процессы кристаллизации (плавления), но и сами испытывают влияние этих процессов. Таким образом, при расчете жидкого моста решается задача Стефана о фазовом переходе (граница жидкость — твердое тело) в области, имеющей свободную поверхность (граница жидкость— газ). Существующие методы численного решения задачи Стефана с конвекцией в областях с фиксированными границами жидкость — газ можно разделить на два типа: методы фиктивных областей (МФО) и основанные на них алгоритмы сквозного счета соответствующих уравнений с разрывными коэффициентами [1-3] и методы с выделением разрыва и точным удовлетворением условий на нем [4]. Совместные расчеты, проведенные авторами [3] и [4] для одних и тех же параметров, показали хорошее совпадение результатов не только для искомых функций (температура, вихрь, функция тока), но и для таких интегральных характеристик, как числа Нуссельта. Основные трудности решения задач в областях, имеющих свободную границу жидкость — газ, связаны с тем, что граничные условия задаются на заранее неизвестных поверхностях, которые должны быть определены в процессе решения задачи. Кроме того, формулировка рассматриваемых при

* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (грант №96-01-01546) и программы интеграционных проектов СО РАН (код проекта №36).

© А. С. Овчарова, 1999.

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

1. Математическая модель I

Пусть вязкая несжимаемая жидкость плотностью р, с кинематической вязкостью V и коэффициентом поверхностного натяжения а(Т) заполняет область 0Б{0 < х < Ь; 0 < у < f (х)}, где х = 0, х = Ь, у = 0 — твердые непроницаемые стенки, у = f (х) — свободная поверхность. Движение жидкости и теплообмен будем описывать системой уравнений Обербека — Буссинеска, которая в переменных ф, ш, 9 имеет вид

дш + д ^ _ д ^ = _1Дш + . (ц)

дЬ дх \ ду) ду \ дх ) Ие Ие2 дх'

Дф + ш = 0; (1.2)

д9 + А (9дф) _ д (9дф) = _^Д9. (1.3)

дЬ дх \ ду) ду \ дх) ИеРг ' '

Здесь Ие, Ог, Рг — числа Рейнольдса, Грасгофа и Прандтля, в качестве масштаба давления выбрана величина ру;2, — характерный масштаб скорости. 9 = (Т _ Т0)/ДТ, где Т0 — характерное значение температуры, ДТ — перепад температур. Вектор силы тяжести g параллелен оси у и направлен вниз.

Функция тока введена соотношениями

дф дф и = —, V = _ —. ду дх

Будем считать, что граничные условия, необходимые для решения системы (1.1) — (1.3), на боковых и нижней границах области ОБ заданы. На свободной границе ставятся кинематическое и динамическое условия. Кинематическое условие имеет вид

Л = V _ и,и. (1.4)

Из динамического условия можно получить в явном виде граничные условия для ш и ф на свободной поверхности. Определим векторы нормали и касательной к свободной

поверхности f (ж) в каждой ее точке как

- Ух 1 I I 1 /ж

п

fx2 л/1 + /2 К/1 + /2 л/1 + /2

и предположим, что коэффициент поверхностного натяжения а(Т) зависит от температуры линейно

1

а(Т) = сто(1 - (т(Т - То)), ( = о(Т0), от =--- |т=То, (,(т > 0.

а0 аъ

Далее введем обозначение. Пусть

8ф дп

уя =

(1.5)

У = / (х)

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

, дvs дР 1 , /х а (л йч

57+^ = - 57 - ке дп + Тгг/!(Оге - с). (1.6)

Здесь и далее время Ъ рассматривается как параметр. О = g0h0/v0 — число Галилея, д0 — ускорение свободного падения, ^ — характерная глубина жидкости.

Условие непрерывности нормальной компоненты вектора напряжений, следуя [5], можно представить в виде

Р - Р = -С" 1 (1 - атАТЛ - -2^ 1= /хх (17)

0 Ке Я V (0 ) Ке бе' Я ^(1 + /х2)3' 1 ' }

Здесь Р — давление жидкости на свободной поверхности; Р0 — внешнее давление (например, атмосферное); Я — радиус кривизны /(ж); Са = /а0 — капиллярное число. Дифференцируя (1.7) по в и подставляя результат в (1.6), получим

. д^ 2 5Ч

+ ^ 57 = ке 5^ + (1.8)

где

1 5ш 5Р0 Са"1 5

— -—--—- - —г- +

Ке 8п 8в Ке 8в

я г! Л - атАТе

Я V (0

+~т(&е - О)-

Из условия непрерывности касательной компоненты вектора напряжений, используя прием, описанный в [5], можно получить явное выражение для вихря на свободной поверхности:

2 мп ее , ,

ш = Я1' + ке 5? (1-9)

Мп = атAT/p0v0v — число Марангони.

Граничное условие для функции тока следует из (1.5), где V' — решение уравнения (1.8). Уравнение (1.8) имеет дивергентный вид. Особенности поведения границы в зависимости

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

Л + у/1+Л дфф = 0. (1.10)

Метод решения. При такой развязке граничных условий на свободной поверхности для решения задачи (1.1) — (1.3) можно использовать любые методы, применяемые для решения задач тепло- и массопереноса в замкнутых областях в переменных ф, ш, 9. Отме-тим,что в рамках модели Обербека — Буссинеска уравнения (1.1)-(1.3) одного типа. Это утверждение относится также и к уравнению (1.2), которое далее будет решаться методом итераций и вводом итерационного параметра может быть подведено под общий тип. Здесь предлагается метод решения в регулярной области. С этой целью отобразим область ОБ на прямоугольник со сторонами 0 < £ < Ь, 0 < п < 1 с помощью преобразования

х = £ у = ] (£ )п.

Тогда все границы области ОБ, включая свободную поверхность, будут лежать на координатных линиях, а каждое из уравнений (1.1)-(1.3) может быть представлено в виде

дФ 1

Ж = в] и (Ф) + (Ф)] + Г

(1.11)

где

Для Ф = ш:

Для Ф = Ф:

дФ дФ дф и(Ф) = БпдФ + В12дФ _ АФдф, д£ дп дп

дФ дФ дф

V(Ф) = В12 * + В22 дП + АФ^

Бц = ] (£), Б12 = _] п, В22 = (1 + В\2)/] (£).

Ог

д9

д9\

Б = Ие, А = Ие, Г = ^ + В12 ^

Б = 1/Л, А = 0, Г = Хш,

(1.12)

Л — итерационный параметр, который вводится при решении уравнения Пуассона для ф. Для Ф = 9:

Б = Рг, А = Рг, Г = 0.

Граничные условия на свободной поверхности (п = 1) примут вид

1 дч.

дvs дvs 2 д /

~ж+Ьа = ие уггт2 д£

+в,

(1.13)

В = 1 (вдш + В дЛ др0 + Са-1 д

В = _ие щ + Б22дп) _+ %

1 Л _ ^9 я V ^0

+ ](Ог9 _ О);

дф = ч

дп

iv/г+]f

Б

(1.14)

22

2 мп ее

ш =о v' + —I-(1.15)

Я к^ 1 + /

/ +1 = 0. (1.16)

Здесь стоит отметить, что кинематическое условие (1.16) теперь выражает закон сохранения массы, что автоматически решает проблему с сохранением объема.

Решение уравнения (1.11) на каждом временном шаге будем искать по схеме стабилизирующей поправки [6], взятой в форме

фк+1/2 _ фк 1

ф т ф = В/ик+1/2(Ф)+кк(ф)]+^

(1.17)

фк+1 - фк+1/2 1

—-— = -щ [кк+1(ф) - ук (Ф).

Здесь

ик+1/2(ф) = ВиФк+1/2 + Б^ - АФк+1/2Зф, ? /

Кк+1(Ф) = -12 фк+1/2 + В22Фк+1 + АФк+1 Ц.

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

Дифференциальные выражения типа (а11ф^, (а22Фч)ч, (а1Ф)^, (а2Ф)ч аппроксимируем со вторым порядком точности конечно-разностными аналогами Л11, Л22, Л1, Л2 соответственно, которые имеют традиционное представление [6, 7]. Для аппроксимации смешанной производной, например (а12Ф^)п, использован оператор Л12Ф = Л2Л1Ф. Оператор Л21Ф определяется аналогично. Тогда схема (1.17) будет аппроксимировать (1.11) с точностью 0(т + После подстановки в (1.17) вместо Ц-(Ф) и Уп(Ф) их разностных аналогов на каждом полушаге по времени для всех внутренних точек получим систему линейных разностных уравнений относительно функции Ф(£п, Пт). Система имеет трехдиагональную структуру с преобладанием диагональных элементов матрицы и может быть эффективно решена методом прогонки с учетом специфики граничных условий. Уравнение (1.13) относится к классу уравнений Бюргерса — Хопфа с правой частью. Второе слагаемое в левой части аппроксимируем следующим образом:

дV' 1 а 2 1- (Ж 1 ып+1ый1 - ып-мП_1 = й"^^') 4 + , + +

—►--1-—► —

2 ее 2 2А£ 2 2А£

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

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

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

В начальный момент времени жидкость, заполняющая область ОБ(0 < х < 1, 0 < у < 1), находится в состоянии покоя. Здесь х = 0, х = 1, у = 0 — твердые непроницаемые стенки; у =1 — свободная поверхность. Начальные условия для температуры заданы в виде

0 < у < 1: 9(х, у) = 0.5 _ х. Граничные условия на твердых стенках:

х = 0: ф = 0, ^ = 0, 9 = 0.5, дп

х =1: ф = 0, ^ = 0, 9 = _0.5, дп

условие Тома. соотношениями (1.14) Vs(1) = 0, Ро = 0. Для

температуры на свободной поверхности задано условие д9/дп = 0. Ниже представлены результаты расчетов, проведенных на сетке 21 х 21, иллюстрирующих эффекты термокапиллярности (конвекция Марангони) и плавучести (конвекция Грасгофа).

На рисунках 1 и 2 изображены изолинии функции тока (а) и изотермы (б) для течения, рассчитанного при значениях параметров

Ие = 1, Ог = 0, О = 0, Са-1 = 1, Рг = 0, 73, Мп =1 и Мп = 2

соответственно; штриховой линией обозначено положение свободной поверхности в начальный момент времени.

На рисунках 3 и 4 представлены результаты расчетов, проведенных для параметров

Ие = 1, Ог = 1200, О = 100, Са-1 = 100, Рг = 0, 73, Мп = 0 и Мп = 40 соответственно.

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

зЗамечание 1. Для реальных физических расплавов величина отДТ9/о0 много меньше единицы, и в условии (1.7), отражающем непрерывность нормальной компоненты вектора напряжений, ею обычно пренебрегают (см., например, [8]). Во втором условии, отражающем непрерывность касательной компоненты вектора напряжений, она сохраняется, так как там она играет важную роль в формировании свободной поверхности. Для корректного сравнения с результатами [8] в приведенных выше расчетах эта величина так же, как и в [8], опускалось.

у = 0 : ф = 0, = 0, 9 = 0.5 _ х. дп

Для связи функции тока и вихря на твердых стенках используется На свободной поверхности граничные условия для ф и ш заданы и (1.15) соответственно, где vs — решение уравнения (1.13); ч.,(0) =

Рис. 1.

Рис. 2.

Рис. З.

Рис. 4.

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

зЗамечание 3. При достижении стационарного решения устанавливался угол подхода свободной поверхности к твердой стенке. Во всех расчетах тангенс угла смачивания (|В12(0,1)| , |В12(1,1)|) не превышает величины 0.02. То есть угол подхода не превышает величины 2°.

2. Математическая модель 2

Пусть СВ = {0 < ж < Ь; 0 < у < f(ж)} — область решения задачи Стефана о фазовом переходе и в* — температура фазового перехода. Граница раздела фаз Б(¿) = {(ж, у) : (ж, у) € СВ, в = в*} разделяет область СВ на две подобласти СВ1 и СВ2, занятые жидкой и твердой фазой вещества соответственно, причем СВ1 = {(ж,у) : (ж,у) € СВ, в > в*}, СВ2 = {(ж,у) : (ж,у) € СВ, в < в*}. f (ж) — верхняя граница области СВ, где крестиками отмечены точки, принадлежащие жидкой фазе вещества, кружочками — твердой фазе (рис. 5). Будем считать, что граничные условия на боковых (ж = 0,ж = Ь) и нижней (у = 0) стенках области СВ заданы; процессы тепло- и массопереноса в жидкой фазе вещества описываются системой уравнений (1.1)-(1.3). Кроме того, для простоты изложения метода предположим, что теплофизические свойства вещества в каждой фазе постоянны; скачком плотности при фазовом переходе можно пренебречь.

В соответствии с методом, изложенным в [3], продолжим задачу (1.1) — (1.3) из области СВ1 в область СВ2 не меняя обозначений:

ди 1 Л — = —Ди, & Ив '

Дф - 1 ф = 0,

£

д

ИвРг

Дв,

(2.1)

(2.2) (2.3)

где 0 < £ < 1 примут вид

малый параметр метода фиктивных областей. Условия согласования

дф дп

Б(1)

0, [и]

в(1)

ди дп

Б(1)

в(1)

условие согласования для температуры есть не что иное, как условие Стефана

'дп

— кт

'дп

И,

которое при сквозном счете уравнения теплопроводности выполняется автоматически.

При таком продолжении функции тока ф сеточные значения этой функции в СВ2 при проведении численных расчетов имеют порядок 10-10 — 10-18, что делает продолжение (2.1), (2.3) вполне оправданным.

Для решения этой задачи использован алгоритм, описанный выше. Как и в ММ1, отобразим область СВ на прямоугольник со сторонами 0 < £ < Ь; 0 < п < 1, а каждое из

1

0

к

У А

О L х

Рис. 5.

уравнений (1.1)-(1.3), (2.1)-(2.3) представим в виде (1.11). Граничные условия на п = 1 примут вид (1.14) - (1.15) для жидкой фазы вещества и ф = ш = 0 — для твердой фазы. Для определения самой свободной поверхности служит уравнение (1.16).

При сквозном счете уравнения (1.11) по методу, изложенному в [3], для всех функций в,ш,ф в области GB нет различия между жидким и твердым состояниями вещества. Граница фазового перехода S(t) определяется автоматически, по значениям температуры в = в* .В случае, когда жидкая фаза вещества имеет свободную деформируемую поверхность (жидкость — газ), ситуация усложняется еще тем, что точка контакта твердой и жидкой фаз, так же как и свободная поверхность, заранее не известна, а вычисляется в процессе решения задачи. При изменении границы фазового перехода часть свободной границы, принадлежащая жидкой фазе, может либо поглощаться твердой фазой вещества (при кристаллизации), либо увеличиваться за счет увеличения общего объема вещества в жидкой фазе (при плавлении) в результате изменения температурного поля. При этом в постановке задачи отсутствует механизм, который регулировал бы процесс подхода свободной поверхности к точке контакта с твердой фазой. Одним из способов для устранения возникшей неопределенности может служить задание угла подхода свободной границы к точке контакта с твердой фазой вещества или (как его еще называют) угла кристаллизации. Поскольку расчет свободной границы (жидкость — газ) вместе с определением границы фазового перехода — наиболее важный момент в решении общей задачи, опишем его более подробно.

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

Этап I. Из уравнения (1.16) определим новое положение свободной поверхности, соответствующее моменту tfc+i = (к + 1)т, отдельно для каждой из фаз: для жидкой фазы так, как это делалось в ММ1; для твердой фазы f = const, так как ф = 0 на линии п = 1. При этом функция f (x) считается всюду непрерывной. По формулам (1.12) находим матрицу коэффициентов B11, B12, B22.

Угол кристаллизации удобнее всего (и, по-видимому, наиболее правильно) задать именно на этом этапе при определении матрицы коэффициентов B11,B12,B22. Если нет иной информации (например, опытных данных), то можно воспользоваться аналогом задачи, рассмотренной в ММ1, и задать угол подхода свободной поверхности к твердой фазе в точке контакта, например, равным нулю. Это можно сделать, полагая равным нулю ко-

эффициент Б12 на линии ц = 1 в точке и1 (см. рис. 5) — последней точке разностной сетки, принадлежащей жидкой фазе. Все остальные коэффициенты матрицы Б^ для всей области ОБ вычисляются по формулам (1.12). Решая уравнение (1.13), находим граничные условия для ф и ш на свободной поверхности для жидкой фазы по формулам (1.14) и (1.15) соответственно. Для твердой фазы, как уже упоминалось выше, ф = ш = 0.

После того, как будет определена матрица коэффициентов Б11,Б12,Б22, для следующего этапа вводится понятие глобальных итераций.

ЭтАП II. На каждой из глобальных итераций уравнение (1.11) решается последовательно для температуры, вихря и функции тока во всей области ОБ сквозным счетом по схеме стабилизирующей поправки, взятой в форме (1.17):

1) определяется температурное поле, где функция тока ф берется с предыдущей глобальной итерации;

2) решается уравнение для вихря ш со значениями температуры, взятыми с предыдущего временного шага и значениями ф с предыдущей глобальной итерации;

3) определяется значение функции тока.

После этого пункты 1-3 повторяются до тех пор, пока не будет достигнута заданная точность для всех функций 9,ш,ф.

На этом переход на новый временной слой завершен. Далее процесс повторяется, пока не будет получено стационарное решение, которое считается найденным, если выполняется условие

тах - С | < еь

где К — заданное число шагов, е1 — заданная точность.

Результаты расчетов. Пусть в начальный момент времени подобласть ОБ! {0 < х < х0; 0 < у < 1} заполнена жидкой фазой вещества, подобласть ОБ2{х0 < х < 2;0 < у < 1} — твердой фазой , хо — граница раздела фаз (на рис. 6, 7 она обозначена вертикальной штриховой линией). Жидкость, заполняющая подобласть ОБ1, находится в состоянии покоя. Начальные условия для температуры заданы по линейному закону в жидкой фазе 0 < у < 1: 9 = (9* — 9\)х/х0 + 9\,

в твердой фазе 0 < у < 1 : 9 = (92 — 9*)(х — хо)/(2 — хо) + 9*, 9Х > 92, 9* — температура фазового перехода.

Рассматриваются граничные условия

х

х у

Для связи функции тока и вихря на твердой стенке х = 0 использовано условие Тома. На поверхности /(х) граничные условия задаются в зависимости от того, к какой фазе вещества (жидкой или твердой) принадлежит эта граница (см. ММ1). Для температуры на границе / (х) задается условие д9/ди = 0. Угол кристаллизации полагался равным нулю. Все расчеты проводились на сетке 41 х 21.

На рис. 6 изображен процесс плавления (а — изолинии функции тока, б — изотермы). Расчет проведен для параметров

0 : ф = 0, дф ди = 0, 9 = 91,

= 2 : ф = 0, ш= 0, 9= --92,

= 0 : ф = 0, ш= 0, д9 ди = 0.

И,е = 1, Ог = 1200, О = 100, Са-1 = 100, Рг = 1, Мп =10,

а б

1.13 1.13

О ж0 2 о х0 2

Рис. 6.

.10 б МП

О ж0 2

Рис. 7.

Хо = 0, 7; вг = 1; в2 = -1, 65.

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

Ив = 1, Сг = 1200, С = 100, Са-1 = 100, Рг = 1, Мп = 10,

Хо = 1; вг = 1; в2 = -2, 25.

Замечание 4- Из рис. 5 видно, что в алгоритм уже "заложено" нарушение закона сохранения объема. При достижении стационарного решения относительная погрешность составляет 0,2% при расчете плавления (см. рис. 6) и 0,025% при расчете кристаллизации (см. рис. 7). Погрешность может быть уменьшена при проведении расчетов на более подробной сетке. Кроме того, можно воспользоваться приемом, позволяющим проводить расчет свободной поверхности на подробной сетке, а расчет внутри области на редкой.

Автор выражает благодарность Ж. Л. Коробицыной и В. В. Кузнецову за полезные консультации и заинтересованное обсуждение результатов работы.

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

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

О 2

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

[3] КОРОБИЦЫНА Ж. Л., Тычков С. А. Численное моделирование процессов тепло- и массопереноса с учетом фазового перехода в геодинамике. ЖВМиМФ, 37, №6, 1997, 733-741.

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

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

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

[7] САМАРСКИЙ А. А. Введение в теорию разностных схем. Наука, М., 1971.

[8] Chippada S., Jue T.C., Ramaswamy D. Finite element simulation of combined buoyancy and thermocapillary driven convection in open cavities. Intern. J. Numer. Meth. Eng., 38, 1995, 335-351.

Поступила в редакцию 1 июля 1998 г., в переработанном виде 8 октября 1998 г.