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

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

Том 10, № 3, 2005

РАСЧЕТ ДВИЖЕНИЯ ЖИДКОСТИ С ПЕРЕМЕННОЙ ВЯЗКОСТЬЮ В ОБЛАСТИ С КРИВОЛИНЕЙНОЙ ГРАНИЦЕЙ

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

Новосибирск, Россия

In this paper an axisymmetric non-stationary flow of a liquid under the combined action of thermo-capillary, ponderomotive, buoyancy forces and rotation is considered in application to the problem of a floating-zone melting in the magnetic field. The conservative monotonous difference scheme and the orthogonal coordinate system is used. A calculation of the convection in the floating zone with the specified shape of its boundary is performed.

Введение

Рассматривается нестационарное осесимметричное движение жидкости, вязкость которой зависит от температуры. Сечение области течения плоскостью ip = const, где ip — полярный угол, имеет форму криволинейного четырехугольника. Предполагается, что известно конформное отображение прямоугольника на этот криволинейный четырехугольник. (В данной работе рассчитывается пример, в котором такое отображение задано аналитически. В общем случае следует использовать численный алгоритм, описанный в работах [1, 2].) Задача ставится в переменных вихрь — функция тока в соответствии с приближением Буссинеска. Дифференциальные уравнения для вихря, функции тока, азимутальной компоненты скорости и температуры записываются в дивергентной форме. Дифференциальные операторы конвективного и диффузионного переноса аппроксимируются монотонными консервативными разностными операторами, имеющими второй порядок аппроксимации по пространству при малых числах Рейнольдса и Пекле и первый — при больших числах Рейнольдса и Пекле (здесь используются результаты работы [3, с. 280-288]). При решении разностных уравнений используется метод, позволяющий точным образом разделить задачи вычисления вихря и функции тока, описанный в работах [4, 5]. При решении на каждом временном шаге уравнения для функции тока используется эффективный итерационный алгоритм, не требующий информации о свойствах и границах спектра разностного оператора, обобщающий алгоритм, описанный в работе [6].

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

© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2005.

малыми скоростями и вращаются в противоположных направлениях. Между ними находится плавающая зона, поддерживаемая в жидком состоянии неподвижным источником высокочастотного электромагнитного поля — индуктором. Токи, наводимые индуктором, сосредоточены в тонком скин-слое, примыкающем к свободной границе расплава. Они создают пондеромоторную силу, направленную ортогонально свободной границе, экспоненциально убывающую при удалении от нее и являющуюся одним из источников конвекции в расплаве. Форма функции пондеромоторной силы на свободной границе задается в данной работе аналитически. Но размеры области и энергетические характеристики процесса, принятые при расчете, соответствуют реальным. На основе разработанного метода решается гидродинамическая часть задачи — производится расчет конвекции в плавающей зоне при заданной форме ее границы.

1. Размерные параметры задачи

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

V(Т) = аи + К/(Т - си).

Ниже приведены названия, обозначения и значения материальных констант [7, 8]. Символ т означает жидкую, в — твердую фазу кремния. Коэффициент аи вычисляется по формуле

ау Vm Ьу/ (Тт Оу )•

Кроме перечисленных нам понадобятся следующие размерные параметры: толщина скин-слоя [9] £т = с/\/2^7тш0 = 2.85 • 10-2см, радиус монокристалла гс = 5 см, скорость протягивания монокристалла ус = -5.52 • 10-3 см/с, угловая скорость вращения верхней Uf = -2.16 рад/с и нижней = 5.24 • 10-1 рад/с границ фазового перехода (соответствующих заготовке и монокристаллу) [7], характерный размер в плавающей зоне I = 2 см,

Название параметра

Символ

Значение

Круговая частота тока Постоянная Больцмана Коэффициент серости Температура плавления Значения р при Т = Тт

Значение V при Т = Тт -(!р/(1Т при Т > Тт -<1а/<1Т

Коэффициенты уравнения для V

Ускорение свободного падения Удельная теплоемкость Коэффициент теплопроводности Скорость света

Коэффициент электропроводности

ио <1

/ьт Т

Тт рт

Рв

Vm

РТ <Т

к

Су

9

Срт ^т

С

7т

107 рад/с

10-5 эрг/(с- см2-Х4)

1.76 5.67 0.27 1700 X 2.53 г/см3 2.30 г/см3 3.2 • 10-3 см2/с 1.52 • 10-4 г/ (см3 •К) 0.1 дин/(см-К) 4, 776 • 10-2 К-см2/е 1661 К 980 см/с2

7, 45 • 106 см2/(К-е2) 5.2 • 106 г-см/(К- с3) 3 • 1010 см/с 1016 1 /с

характерная скорость движения расплава У\ = 10 см/с [10], характерный перегрев расплава ДТ = 50 К [7], максимум напряженности магнитного поля на свободной границе Н0 = 300 г1/2/(е-см1/2) [10]. Будем считать, что радиус заготовки г/ = гс. Тогда и скорость протягивания заготовки V/ = vc.

2. Безразмерные критерии подобия

Выберем в качестве масштабов длины, скорости, времени, вязкости и температуры соответственно l, vi, l/vi, vm и AT. Тогда в задаче появятся следующие безразмерные параметры:

Re = ^ = 6250, Gr = = 2.300 ■ 106,

Vm pmVm

H02l ат ATI 5

E«m =-0-2 = 993.3, Ma = ——— = 3.860 ■ 10

8n£mPrnV2 Pm^

Vc = — = -5.520 ■ 10-4, Vf = f = -5.520 ■ 10-4,

v1v1

f = f = -0.4320, = ^ = 0.1048, S = = 0.9C

v1 v1 Pm

Pe = = 72.49, Qe =-= 13.37,

Am 8npmCpmViAT

Bi = /A^Ar = 0.9836, To = aAt = 34,

a b c

Av = — = 0.6173, Bv = —^ = 0.2985, Cv = = 33.22, Vm ATvm AT

Яс = у = 2.5, Я/ = г/ = 2.5, Ет = у = 1.425 • 10-2.

Здесь Ке — число Рейнольдса; Ог — число Грасгофа; Еит — магнитное число Эйлера; Ма — число Марангони; V:, V/ — безразмерные скорости протягивания; Пс, П/ — безразмерные скорости вращения; Б — отношение плотностей; Ре — число Пекле; Qe — отношение характерных интенсивностей джоулева тепловыделения и конвективного теплопереноса; Б1 — число Био; Т0 — безразмерная температура плавления; Аи, Ви, Си — коэффициенты уравнения для безразмерной вязкости

V = Ау + В/(Т - С);

Яс, Я/ — безразмерные радиусы монокристалла и заготовки; Ет — безразмерная толщина скин-слоя.

3. Уравнения и граничные условия

Пусть Or, z, р — цилиндрическая система координат, где r — полярный радиус, р — полярный угол. Рассмотрим сечение р = const. Направим ось Or вправо, а ось Oz — вверх. Пусть функции r(x,y), z(x,y) осуществляют конформное отображение прямоугольника 0 < x < 1, 0 < y < Y на область, занятую расплавом, так, что стороны x = 0, x = 1,

у = 0, у = У прямоугольника переходят соответственно в ось симметрии (г = 0) Г0, свободную (правую) границу Гт, фронт кристаллизации Гс (нижнюю границу расплав — монокристалл) и фронт плавления ^ (верхнюю границу расплав — заготовка).

Пусть п(х,у,Ь), у(х,у,Ь), /ш(х,у,Ь) — компоненты скорости расплава в направлениях х, у, <£>. Введем функцию тока Ф и вихрь ш по формулам

u

1 д Ф

rH ду '

v = —

1 д Ф

rH дх '

и

1

H2

д(иН) д^Н)

ду

дх

(1)

где

Н = \J (дг/дх)2 + (дг/ду)2

— коэффициент Ламэ.

Форма уравнений осесимметричного движения жидкости с переменной вязкостью следует из тождеств

11 — (divP)w =

pm

д ( 3 д fw\\ д f 3 д /w V r2H2 \, дх V дх \r

ду ду r

где

1

r2H 2

д ( дт- \ д ( дт--— 2v—rw\ - — 2v—rw дх дх ду ду

+

д д д д + дХ KvYx(rw)) + ду У^дТу(rw)

1

PmH 2

д^ (divP)x) д^ (divP)y)

дУ

дх

1

H2 1

H2

дх

д дх

1íruu)\ + А (1íruu)\ + дА + дВ

r дх J ду \r ду J дх ду

2 дr ди\ и\ д

дх дх r ду

2 дr ди\ и ду ду r

д ( д (и \\ д ( д (и \\ дА дВ

+--rv— — +--rv— — +---1--

дх дх r ду ду r дх ду

= д^ 2 ди 2 дH дх \ H ду H2 дх

= д V (ди А дЕ

дх \H дх H2 ду

дv Í 2 дv 2 дH

ду \H ду H2 дх

дv ( 2 дv 2 дH

ду \ H дх H2 ду

u

u

+

P — тензор напряжений. Эти тождества устанавливаются прямым вычислением. При этом используются формулы для компонент P согласно [11] и компонент divP согласно [12]. Формулы проверены в пакете "Математика" для конкретного случая, когда

x x • т т x

r = e cos у, z = e sin у, H = e .

Введем модифицированные скорости и, V, Ш и модифицированный вихрь П по формулам

Ш = гт, П = Ш, (2)

г

и = дФ, V = - дФ.

ду дх

(3)

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

где

дП 1 Ж + гН2

д (( 1 Л дг д/. и I п )

дх V V Ке V дх дх

+ д / / 1 / дг+гдv^+v\ п) ду \\ Ке V дУ ду/

(дП

Ке \дх \ дх у ду \ ду

С1

С1 =

гН2

д /дгЖ2

дх уду г3

д

+

дг W2

+

ду \дх г3

1 дНд Ф

(1 (дФ) _

Ке \ дх \ дх \ Н ду \ гН ду / гН3 дх дх

д^ / 1 д / 1 дФ\ 1 дНдФ ду \ Нду \гН дх у гН3 дх ду

1 д Ф\ 1 дНд Ф

- -

ду \дх \Ндх \гН ду у гН3 ду дх

д^ / 1 д ду

Ог / д

1 д Ф\ 1 дНд Ф Н дх V гН дх у гН3 ду ду

+

Бе2 V дх V дхТ) + ду (дуТ) ) Е^тду (Н/„)_

у

/„(х,у) = /„I [ Н(1,у)^у/вМ ехр(- 2Н(1,у)(1 - х)/Ет)

(4)

(5)

безразмерная пондеромоторная сила [9];

г

в = Н (1,у)^у;

/„(х) — функция пондеромоторной силы на границе Гт, имеющая единичный максимум и заданная на единичном отрезке.

Выведем граничные условия для П. На оси симметрии должно выполняться условие симметрии

дП

(6)

дП

— = 0, х = 0. дх

На свободной границе Гт — условие Марангони [13, с. 195]

- ™ - да

в • Р • п = —, дв

1

где п — внешняя нормаль к границе; в — вектор, получающийся поворотом п на 90° по часовой стрелке, или

_ 1 да _ ат дТ

-Рху _ -ИдУ _ Иду- (7)

Согласно [11]

/1 ди 1 ду и дИ V дИ\

Рху _ р^улдУ + ИдХ - И2 ~дУ - н2 ~дХ) - (8)

Используя (1), (2) и условия и _ ди/ду _ 0 на Гт, представим рху в виде

/ ^ 2у дИ \

Рху _ -рт и + — ~дх),Х _ 1-

Подставляя это выражение в (7), обезразмеривая и пренебрегая членом, содержащим множитель у, ввиду большого отношения параметров Ма и Ке, окончательно получим

П _ Ма 1 дТ _ 1 Ке г Ни ду '

Член, содержащий множитель у, мы отбросили для того, чтобы в дальнейшем точно расщепить задачи вычисления для П и Ф.

На границах фазового перехода поставим условия проскальзывания [14, с. 129]

в■ Р ■ п _ рху _ кв ■ (V — и), у _ 0, у _ У,

где V — скорость стенки (твердого кремния); и — скорость расплава; к — коэффициент проскальзывания. Используя (1), (2), (8), представим рху в виде

/ 2 ду 2идН

рху _ + НдХ — И -дУ

Компонента у скорости и на Гс,/ равна проекции скорости стенки на направление у, умноженной на Б:

_ 1 дг у _ ус,/э — ^. И ду

Используя определение Н и условия Коши — Римана, находим, что

_ —Кс/ ,

ду ^дг 1 дН

— _ —ус,/^Кс/

дх дх ' Н2 ду

где

- с,/

(10)

К 1 / д2г дг д2г дг

с/ Н3 у дх2 дх дх2 дх

— кривизна границ фазового перехода. Кроме того,

^ ^ . ус/ дг _ _ в ■ у _ %с ^^^—, в ■ и _ гс /и,

с/ Н дх с/

где гс _ —1, г/ _ 1. Проведя обезразмеривание, получим условия проскальзывания в следующем виде:

1 дФ V дг

П _ -(2К — А)^ +(2КЭ — АА17НШ' у _ 0

Q

1 дФ Vf

-(2K + тщ+ <2KS+Ai) тЯ d;-y = Y

(ii)

где Kc,f определено в (10); А/ — безразмерный параметр проскальзывания,

а = -kL.

pmVm

(В дальнейшем при расчетах полагается А/ = 100.) Здесь учтено, что на Гс,/ безразмерная вязкость V =1.

Отметим, что постановка условий проскальзывания вместо условий типа прилипания обусловлена тем, что последние вызывают разрыв вектора скорости в точках трехфазного контакта, так как протяжка осуществляется параллельно оси г, а свободная граница в этих точках не параллельна оси г. При численном решении задачи такого разрыва не возникает, так как вместо точных используются приближенные условия типа условий Тома. Можно показать, что главные при к ^ 0 члены в этих условиях определяют условия проскальзывания с непостоянным А/ = 1/(кН), где к — шаг сетки в направлении, ортогональном к рассматриваемой границе.

Задача для функции тока

(\дФ\ + дФ\ = гЯ2п

дх \т дху ду \т ду у '

(12)

Ф

0, х 2

0, Ф

R2

--f VfS, х =1,

Ф = - у VcS, у = 0, Ф = - у Vf S, у = Y. Функция W должна удовлетворять уравнению импульса

(13)

дW 1

+

öt тЯ2

+ А ИЛ + Ч + VI WI -

дх \ \ Re дх у у ду \ \ Re ду

1

д

, , vt-

Re дх дх

öW

д / öW\ ду V ду У

0,

граничным условиям

W = 0, х = 0

(14)

(15)

(это следствие из (2)), условию отсутствия касательного напряжения в направлении р на свободной границе [11]

1 д(W/т) W дт

Pxf

или

Я дх

öw 2 ÖTW

дх дх r

r2H дх 0, х = 1

0

и условиям прилипания на твердых границах

W = Псг2, у = 0, W = П/г2 Для Т должны выполняться уравнение энергии

Y.

дТ 1

Ж + ТЯ2

д(иТ) + ö(VT )

дх

ду

1

дх

дТ д дТ

т— +--т—

дх ду ду

Qefn

(16)

(17)

(18)

у

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

дТ

— _0, х _0, (19)

дх

закон Стефана — Больцмана на свободной границе

1 дТ „ / Т 4 4

нэх+вч ть; =0-х=1 (20)

и условия первого рода на фронтах плавления и кристаллизации

Т _ То, у _0, у _ У. (21)

Зададим также начальные условия

Ф _ Фь, П _ Пь, Ш _ Шь, Т _ Ть, г _ 0. (22)

Итак, задача (3), (4), (6), (9), (11)-(22) служит для определения неизвестных функций Ф, П, и, V, Ш, Т. Компоненты скорости и, у, восстанавливаются с помощью первых двух равенств (1) и первого равенства (2).

4. Численный алгоритм решения задачи

Введем обозначения

дг

ди\

Р1 _ П, р1 _ -— 2—и + г— , д1 _ - — 2—и + г—

Ке дх дх

дг

ди\

Ке ду ду

1 _ — Р2 _ Ш 2 _— — и 2 2 _ ^

Ке Ке дх Ке ду Ке

С2 _ 0, Р3 _ Т, р3 _ 0, д3 _ 0, р3 _ С3 _ де/„.

Ре

Тогда уравнения (4), (14), (18) можно записать в универсальной форме

(23)

д

— + Ь\ + Ь12) Г _ С, г _ 1, 2, 3,

(24)

где

Ь\Р *

1

гН 2 1

д((р* + и)Р*) д ( iдР

дх

дх дх

ЦР * _ —-2 гН2

д((д* + V)Р*) д ( IдР*

/I

ду

ду ду

Пусть задана неравномерная сетка

0 _ хо < х1 < . . . < хм-1 <хм _ 1, 0 _ уо <у1 < ... < ум-1 <ум _ У.

(25)

Обозначим

хп-1/2 _ (хп-1 + хп)/2, кхи-1/2 _ (х„ - хга-1), п _ 1, N

ут-1/2 = (ут-1 + ут)/2, кут-1/2 = (ут - ут-1), т = 1, М

кх„ = (х„+1 - х„-1)/2, п = 1, N - 1,

кут = (ут+1 - ут-1 )/2, т = 1, М - 1.

Будем считать, что

кж„+1/2 - кж„-1/2 = 0(кХ„+1/2) , кут+1/2 - кут-1/2 = 0(к'2т+1/2).

Введем также равномерную сетку на полуоси £ > 0 с шагом т. Примем для любой функции /г(х,у,£) обозначение

/ (х„,ут, кт) Угат,

причем индексы п, т могут быть дробными, а индексы п, т, г или к — отсутствовать, если / не зависит от х, у, г или £ или когда они несущественны. В работе [3] для дифференциального оператора

I Р = и — - —

дх дх дх

получено разностное представление

^ и„ -1/2 (Р„ - Т„-1) и„+1/2 (Р„+1 - Тга) г

„ = 2 - 1 2 - „

' ¿'ХП ^ ' ¿'ХП

где

,1 ^„+1/2(Р„+1 - Р„ )

а„+1/2 СШ а„+1/2 —

к

-а„-1/2 еШ а„-1/2 а„-1/2

кж„+1/2

1/2(Р„ - Т„-1) кж„-1/2

2^

Рассмотрим дифференциальное тождество

„—1/2

д(кР) - дР Яи

—--- = к--+ Р—.

дх дх дх

Аналогично ему имеет место разностное тождество

1 (р„+1 + т„) — (р„ + Т„-1) к~ Г„+1/2 2 и„-1/2 2

= и-„+1/2 (Р„+1 - Т„) + ии„-1/2 (Р„ - Т„-1) + р (и„+1/2 - ии„-1/2) = 2 к + 2 к + р„ - '

Исходя из этих двух тождеств для дифференциального оператора

ьр = д (ир) - / дР ^

дх дх дх

(26)

получим следующее разностное представление:

Л Р _ 1 (и (Рп+1 + Рп) и (Рп + Рп-1) \ г ЛРп _ - Vип+1/2-2--ип-1/2-2-у — In,

где 1п определено в (26).

С учетом этого аппроксимируем дифференциальные операторы Ь\, Ьг2 разностными операторами А\, Л2 следующим образом:

Лкр*к 1 \{ттк | „Лк \ р*к

А гк тргк _ ___г (ттк , гк \и*к _

Л1 Рпт _ л2 А [(ип+1/2т + рп+1/2т) Рп+1/2т

гпт Ипт1хп

Рп+1/2тап+1/2т ага+1/2т ^ ^ хп+1/2

_(и к + глк )Р*к ^п+1/2т п+1/2т_п+1/2т ( Р *к _ Р*к ) +

(ип-1/2т + рп-1/2т)Рп-1/2т Н ' (Рп+1т Рпт) +

+ Рпк-1/2тап-1/2т /2т ( Р*к _ Р*к )]

+ 7 (Рпт рn—1m)],

1хп-1/2

п _ 1,N— 1, т _1,М — 1, г _ 1, 2, 3,

Л*к Р*к _ _1_ГП/к + д*к ) Р*к _

Л2 Рпт _ и2 А [( 'пт+1/2 + дпт+1/2)Рпт+1/2

гпт И пт Аут

Рпт+1/2впт+1/2 впт+1/2 у Аут+1/2

_(Vк + д*к ) Р*к <-~пт+1/2 пт+1/2 ' пт+1/2 ( Р*к _ Р*к ) +

пт-1/2 + дпт- 1/2)Рпт-1/2 1 (Р пт+1 Рпт) +

+ Рпкт-1/2впт-1/2 впт-1/2 (Р*к _ Р*к )]

+ 7 (Рпт рnm—1)l,

Пут-1 /2

где

п _ 1,п1, т _ 1,М— 1, г _ 1, 2, 3,

ап+1/2т _ Ахп+1/2(ип+1/2т + рп+1/2т)/(2рп+1/2т);

впт+1/2 _ 1ут+1/2(^„т+1/2 + дпт+1/2)/(2рпт+1/2); (27)

п1 _ N — 1; п2 _ п? _ N.

Распространим операторы Ь\, Ьг2 и Л^, Л2 при г _ 1, 3, когда при х _ 0 ставятся условия второго рода (6), (19), и функцию С1 на случай х _ 0. Переходя к пределу (по правилу Лопиталя) при х ^ 0, в (25) с учетом р _ 0, и _ 0, Рх _ 0, Н _ гх при х _ 0 получим (нижний индекс обозначает взятие производной по соответствующей переменной)

^Р* _ г3 Г((рх + их)Р*)х — ((2рх — р*)Рх)х],

гх3

^Р* _ г3Г((дх + Vx)Р*)у — (РхРУ)у], х _ 0, г _ 1,3.

гх

Продолжим сетку хп на один узел влево и учтем, что функции рх, их — нечетные, а функции рх, р, Р, гх — четные по х. Получим следующую разностную аппроксимацию оператора Ь1:

гкргк 2

\*к р Л1 Р

0т г3

гхот

(рк + ик ) Р*к (2р*к _ р*к )(Р*к _ Р*к )

(рх 1 / 2т + их1/2т)Р1/2т (2рх1/2т р1/2т)(Р1т Р 0т )

Ьх1/2 Ь2х1/2

Аппроксимация оператора Ь2 при х = 0 имеет вид

1fc 7Trtfc

Д ik 77'

Л Fr

Ж'

3 h [(Vx0m+1/2 + qx0m+1/2) ^Om, +1/2

_(Vk I ) pik _ f"x0m+1/^1m+1/2 """ '"1m+1/2 (pifc _ Fife ) +

(Vx0m-1/2 + Уж0т-1/2)p0m-1/2 h - (F0m+1 F0m) +

^ж0т+1/2в1т+1/2 cth в1т+1/2

ym+1/2

1 /2в1т-1 /2 cth в1т-1 /2 hym-1 /2

(pm- F0m-1)], m=I,M- i, i = 1,з.

Функция G1 при X = 0 имеет следующий предел:

1 2 Gr

G1 = — [-/1 + /2 + Re(-/3 - /4 + /5 - /б + /7 - /s - - 1ю) - R2 (I11 + I12)],

где

W2 r 1

г _ Wxx'yx г _ 1

11 = ' /2 =4

W2

T _ - xx

^3 = -

rx

Ф

yxx

2vxxrxxx^xx

v

yxx

Фх

+

rx

1Фх

4 Ф

з rx2

з rx3

'6 =

vy гхххФухх

7=

^жжФужж

vxxrxy Фхх

v / 1 Ф

y I xxxx

4 Фхх rxxx

3 r2

з r3

10

vy rxy Фухх

/11 zy Txx + zyxx/12 [zxxT]y •

Причем с точностью 0(х1)

r1m x1

6r

1m

+

6r

2m

X1 (Х2 X1) X2 (X2 X 2 )

Фх

2Ф

1m

Фх

24Ф

1m

+

24Ф

2m

X

2 2 2 2 2 2

X1(X2 X1) X2(X2 X1)

Wx

vx

2(v1m - v0m)

2(z1m z0m)

2Wi

1m

X

T

x

2(T1m - T0m)

2 2 xx 2

/y»^

1 1 1

Для получения этих формул достаточно разложить все функции, входящие в правую часть (5), в ряды Тейлора по х и перейти к пределу при х ^ 0.

Граничные условия (16), (20) запишем в универсальной форме

dF1 p1 . (i - 2)BiH - - -p — -

SX ^ T04

F14, X = 1, i = 2, 3.

С учетом этого распространим операторы Л1, г = 2, 3, на случай п = N [15, с. 45]:

2

Af F^m = тт2 j [ (UN — 1 /2m + pNv-1/2m)FN-1/2m +

r«mHNmftx«-1/2

4

2

2

4

r

r

r

x

x

x

y

y

vyrxxx

5

2

2

r

r

r

x

x

x

y

y

1

/

8

4

3

4

2

r

r

r

x

x

x

y

y

9

4

r

r

x

x

y

y

r

r

x

, (i - 2)BiHNmßNm piM , PN-1/2maN-1/2m cth ajV-1/2^ ik ^

+ гр4 F Nm + J (FNm F N-1m)],

T 0 hxN-1/2

i = 2, 3.

Аппроксимируем (24) следующим образом:

_ ( FiV+1 _ FiV ) _

/I? I „,2^2 X ik \ik\ (F nm Fnm) , /X ik , AikW77«k+1 pik \ ,

(E + Y T TV 1 Ä2 )- + Y(Л 1 + Ä2 ) (Fnm - Fnm) +

T

ik ik ik ik

+ (Л1к + Л2к)Рпт _ 0гпт, п _ по, п1, т _ 1, М — 1, г _ 1, 2, 3, (29)

где Е — тождественный оператор; Л 1к — оператор, отличающийся от Л1к только при г _ 3, п _ N:

_ ,к 2 к 'к -'к Л1 СМт _ 772 7 (им-1 /2т + рЛк-1/2т)Слк-1/2т+

гМтИМт1хМ-1/2

: 4(г — 2)В1НмтРМктРМк,т Лк , РМу-1/2таМу-1/2т аМУ-1/2т

"ЧМт + 7 Х

1 rp4 Wm 1 7

T0 hxN -1/2

Х (CNm — C_iV-1m)\, i = 3;

Y G [0,1\ — весовой коэффициент; n1 = Щ = 0, n0 = 1.

Уравнение (29) можно переписать в факторизованном виде:

_ ( Fik+1 _ Fik )

/I? I „„Л ik\(p I „„Aik\ (F nm Fnm) /~tik (E + YTЛ1 )(E + YTЛ2 )- = Gnm-

— (Л1к + Л2к)Р*т, п _ п0,п1, т _1,М — 1, г _ 1, 2, 3.

Поэтому переход с к-го на (к + 1)-й временной слой при фиксированном г может быть осуществлен по приведенному алгоритму: 1. Определяем невязку

ппкт _ спт — (Л1к + Л2кп _ п0,п1, т _ 1, М — 1.

2. Для всех т _ 1, М — 1 решаем системы уравнений относительно (_пт:

_0кт _ 0 при г _ 2,

(e+YTÄf )znm = nnm, n = no,n1,

(gik+^ gik ) xik / rp i \ik\\Jim У 1m/ „ • i

ZNm = (E + Yt^2 )- при « = 1

T

где g1m — аппроксимация правой части (9).

3. При всех n = n0, n1 решаем системы уравнений относительно поправки (n^m:

(gik+1 _ gik )

/■ik _ (g2n g2n) /nn\

Sn0 = ; (30)

T

(E + YT&2kKikm = (iL, m = 1,M - 1;

ik (g3n+ - #3^

>ik _ \a3n_a3nJ /01 \

SnM = , (31)

T

где #2„, — аппроксимация правых частей (11) при г = 1, (17) при г = 2 и (21) при г = 3. Далее полагаем

(^ - ^ )

(¿т = 0, если г = 2, <^т = —-—, если г = 1, т = 0, М.

т

4. Для всех п = 0, N т = 0, М вычисляем

ргк+1 _ ргк | гк

р„т р„т + ' Ц„т.

Описанная разностная схема является монотонной при достаточно малом т и консервативной [3, 16]. Если функции и, V и одного порядка, то она превращается в схему с центральными разностями для конвективных членов, которая имеет второй порядок аппроксимации по пространству. Если и, V ^ то она близка к схеме для уравнений Эйлера с разностями против потока и имеет первый порядок аппроксимации по пространству. По времени описанная схема имеет первый порядок аппроксимации.

Обратимся теперь к способу реализации алгоритма нахождения функций П, Ф, при котором разностные аналоги условий (11) выполняются точно. Представим П в виде (несколько иначе, чем в [4, 5])

П„т = П„т + Р„т ($2„+ - $2„) + Qmn ($3„+ - 9ш), (32)

где функции Р^1, ^„т1 удовлетворяют условиям

(Е + 7тЛ2к )Р„т1 = 0, Р„0+1 = 1, РкМ = 0, (Е + ттл2 к)д„т1 = 0, д„+1 = 0, д„М = 1.

Тогда при реализации алгоритма для П следует положить Р„т = П„т, Р„т+1 = 6„т1, а правые части формул (30), (31) при г = 1 заменить нулями. Разностная задача для Ф выглядит следующим образом:

„х фк+1 + фк+1 + ау фк+1 + ьу Фк+1 — г фк+1+

и„т „— 1т "Т" °гат „+1т "Т" и„т „т—1 "Т" °гат „т+1 °„т*„т "Т"

+<т+1Ф„+1 + е„кт+1Ф„М1-1 = -/„т1, п = ^Г-Г, т = 1ГМ-Т; (33)

о2 _

Фот1 = 0, Ф*+1 = -^ т = 0, М;

г2 г2

Ф„+1 = Ф„+1 = V;5, п = бТЖ (34)

где

= -;

„т ь — ' „т г^ — '

г„- 1/2т кх„—1/2 кх„ г„+1 /2т кх„+1/2 кх„

ау = _1_=_; ЬУ = _1_;

„т „т

7 — ' „т 7 — '

г„т-1/2 кут — 1 /2 кут г„т+ 1/2кут+1/2 кут

г = + + аУ + ЬУ ;

„т „т „т „т „т

^ук+1 = _ Р„тге Н„тг„т(А/ - 2Кс„) ; ^ук+1 = _ ^„т Н„тг„т(А/ + 2К/„) ;

„т Н„0г^0ку1/2 ' „т Н„М г^м куМ —1/2

/ = _г Н2 (6^ _ Рк+1^1^ Пк+1«1к) - ^Ук+1Фк+^ рУк+1фк+1_ у„т ' „тН„т(П„т Р „т У2„ ^„т У3„) а„т Ф„0 е„т Ф„М

rnmНПт^сРП)+1(2Ксга$ - A1)(zn+10 - zn-10) 2 hxnrn0

rramH^mV/Qnm1(2K/nS + A1)(zn+1M - zn-1M)

2к х„ г„М Н„М

Для решения задачи (33), (34) используем подход, предложенный В. Г. Зверевым в работе [6]. Сначала рассмотрим так называемый полинейный метод решения данной задачи (для простоты записи опустим верхний индекс к + 1):

Ф+1/2 _ г ф1+1/2 + Ф+1/2 = _ ау Ф+1/2 _ Ьу Ф

""„т „— 1т °„т „т °гат „+1т ""„т „т— 1 °гат „т+1

^„тФ„11/2 - е„тФ„м-1 - /„т, п = ТУЖ-Т, т = 1,М - 1; (35)

аУ ф1+1 _ г ф1+1 + ЬУ Ф1+1 + ^У Ф1+1 + рУ Ф1+1 =

"„т „т-1 °„т * „т 1 "„т „т+1 1 "„т „1 1 ^„т „М-1

= -«„тФ„+-11т - Ь„тФ„++1/,т - /„т, т = 1,М - 1, п = 1,Ж - 1. (36)

Недостатком этого алгоритма является то, что члены Ф„т+1, Ф„м-1 в формулах (35) и член Ф„+1т в формулах (36) берутся с нижнего итерационного слоя, что сильно ухудшает сходимость алгоритма. Поэтому, следуя [6], предположим, что

Ф„т+1 = а„тФ„т + /6„тФ„1 + 6„тФ„М-1 + 6„т; (37)

Ф„+1т '6„тФ„т + п„т. (38)

Подставляя эти выражения в (35), (36) и учитывая неявно Ф„т, Ф„1, получим

ax Ф+1/2 _ (c _ by a W+1/2 + bx Ф+1/2 = _ay Ф+1/2 n-1m V°nm unmLXnmJ nm "T" °ram n+1m ""ram nm-1

(dn,m + bnm/^nm) Фп1 (^m + bnm7nm^nM-1 (/nm + ^m^nm^ (39)

ay Фг+! _ (c _ ьх £ )Ф1+1 + ьУ Ф1+1 , + dy Ф1+1+

anm Фnm-1 (cnm bnm£nmi*nm + bnmTnm+1 + dnmTn1 +

+ermФnM-1 = - anmФn—1m - (/nm + bnm^nm.^ ). (40)

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

anmФnm—1 (cnm 0(anm + bnm))Фnm + bnmФnm+1 + dnmФn1 + enmФnM -1

-anmФn-1m - bnmФn+1m + 0(anm + ^nmi^nm - /nmi (41)

anmФn—1m (cnm 0(anm + bnm))Фnm + bnmФn+1m anmФnm—1

-bnmФnm+1 + ^(anm + ^ram^nm - dnmФn1 - enmФnM-1 - /nm, (42)

где 0 — коэффициент нижней релаксации, 0 < 0 < 1.

Подставляя в (41), (42) выражения (37), (38), получим рекуррентные формулы для определения коэффициентов связи:

anM-1 = ДгМ-1 = ° 7nM-1 = -1 = ФПМ ,

anM-2 = ВДМ - 1J AnM-2 = «M -

1

7nM-2 = 1 + kn(-спм-1 + ®пм-1), -2 = ВДм -am-1 + /)п'м-1),

ау л ьу У + с/у Ьу в + еу

а _ _ пт__У _ иптУпт I 11 пт у _ пт 1пт ' °пт

апт-1 ' 1~.у л , Рпт-1 ' 1~.у л , ¡пт-1 ' ,у л :

Спт Ьптапт Спт Ьптапт Спт Ьптапт

Ьу У + / 'I

У _ пт пт '

пт

пт- 1

Спт Ь1тапт

т _ М — 2.2, п _ 1, N — 1, (43)

где

к-у __1_; /1 _ ах Ф1 , + Ьх Ф1 , - в(ах + Ьх )Ф1 + / ;

п ' ' ^ пт пт п-1т 1 пт п+1т \ пт 1 пт/ пт 1 упт-.

СпМ -1

Спт Спт 0(апт + Ьпт); (44)

У _ П Л+1/2 _ Т У _ аМ-1т

уМ-1т _ 0; Пм-1т _ ФМт; уМ-2т _ С-;

СМ -1т

Ьх П+1/2 + /'1+1/2 ах

Л+1/2 _ ЬМ-1туМ-1т + / N-1т ; У _ _а

ПМ-2т _ II ; уп-1т _

пт

С" ' СII — Ьх У

N -1т Спт Ьптупт

Ьх П+1/2 + /'1+1/2 ___

пп+ут _ ЬптПпт + /пт , п _ ж—22; т _ 1,М — 1. (45)

С'' — Ьх У

пт - пт пт

Здесь

/111+1/2 _ ау Ф1+1/2 + ьу Ф1+1/2 _ в(ау + ьу )Ф1+1/2 +

^ пт пт пт-1 1 "пт пт+1 ^пт 1 пт/ пт 1

+су Ф1+1/2 + еу Ф1+1/2 + / ; С'' _ С _ Л(ау + ьу )

Заметим, что выбор коэффициента кпу в формулах (43) находится в нашем распоряжении. В формулах (44) он определен так, чтобы апм-2, впм-2, Упм-2 и упм-2 вычислялись по тем же формулам, что и апт-1, (Зпт-1, упт-1 и упт-1 при т < М — 1.

Запишем алгоритм перехода с 1-го на (/ + 1)-й итерационный слой при решении системы

(33), (34), I _ 1, 2,...,Ь.

1. Если I _ 1, то по формулам (43) определяются коэффициенты апт, /Зпт, упт, по формулам (45) — коэффициенты Упт.

2. По формулам (43) определяются 5\1т.

3. Для всех т _ 1, М - 1 решаются системы уравнений (39) с учетом граничных условий (34).

4. По формулам (45) определяются Пп+У2.

5. Для всех п _ 1, N — 1 решаются системы уравнений (40) с учетом граничных условий

(34) [17, с. 90].

Описанный алгоритм на порядок эффективнее полинейного метода, так как в нем члены Фпт+1, Фп+1т уже учитываются неявно и только Фпм-1 — явно. Оптимальное значение коэффициента 0 можно подобрать экспериментально.

После решения задачи для Фпт1 вычисляются функции д^^1 и д^1 и по формуле (32) восстанавливается функция Ппт1.

5. Тесты

Если некоторая дифференциальная задача имеет стационарное решение и аналогичным свойством обладает ее дискретный аналог, то можно исследовать дискретное решение на сходимость к решению дифференциальной задачи. Если е1 и е2 — нормы погрешности

какой-либо из искомых функций дискретного решения на сетках N х М и 2Ж х 2М соответственно, то порядок сходимости а для этой функции определяется формулой Рунге

а = ^2(в1/£2)-

Если точное решение дифференциальной задачи известно, то погрешность решения понимается в ее обычном смысле. В противном случае под е* здесь следует понимать норму разности между дискретными решениями на сетках ^ х ¿М и 2^ х 2гМ. Таким образом, в первом случае нужно провести минимум два, а во втором — минимум три расчета (на сетках N х М, 2N х 2М и 4N х 4М). В первом случае можно говорить о сходимости к точному решению, а во втором — о сходимости "в себе".

Рассмотрим следующую задачу:

дФ д П д Ф д П 1 ( д 2П д2П \ ^

+ = ^

ду дх дх ду Ке \ дх2 ду

д2Ф д2Ф

+ тг^ = П при 0 < х < 1, 0 < у < 1, дх2 ду2

дФ

Ф = — = 0, у = 0 и у =1, Ф = 0, х = 0 и х =1, ду

д П

"д~ = 0, х = 0, П = #(у), х =1.

Функции С(х, у) и д(у) подберем так, чтобы функция Ф = Фу = йш2 пх вт2 пу являлась точным решением этой задачи. В дискретном аналоге вместо условий дФ/ду = 0 поставим условия Тома [4]:

22 П(х, 0) = Ф(х,у1), П(х, 1) = —-Ф(х,ум-1).

-у1/2 -уМ -1/2

Эта задача решалась на последовательности сеток методом установления по времени с помощью разностной схемы, описанной в предыдущем разделе. Для каждой из функций Ф, П, и = дФ/ду и V = —дФ/дх вычислялись две нормы погрешности Д (норма в С и в ¿2):

(г 1 У/2

е = тах и а* = I / ^х / Я^у I , г = 1, 2.

\0 0 /

При Ке = 0.1 в обеих нормах имел место второй порядок сходимости, при Ке = 100 — первый. Причем в последнем случае он начинал проявляться при больших числах разбиений, чем в первом. При Ке = 10 000 численное решение не сходилось к точному, так как последнее становилось неустойчивым, но наблюдалась сходимость "в себе" с порядками, лежащими в пределах 0.6 ... 1.0 для всех функций, кроме П. Для П порядок, вычисленный по нормам е*, был отрицательным, а порядок, вычисленный по нормам а*, равен 0.4. Эти результаты получены при 4N = 4М = 160. Отличие порядков сходимости для некоторых функций от 1 при Ке = 10 000 может быть объяснено недостаточным числом разбиений и разрывностью функции П в угловых точках.

Был также произведен расчет порядков сходимости "в себе" а£ и аа, соответствующих нормам е и а* погрешности в определении функций Ф, П, Ж, Т, и, V при решении задачи,

Таблица!

х 4М Параметр Функция

ф п W Т и V

80 х 80 а£ 1.76 1.94 1.95 1.82 1.65 2.02

аа 1.93 1.91 2.03 2.02 1.76 1.92

80 х 160 а£ 1.90 2.32 1.87 1.99 2.01 1.50

аа 1.95 1.94 2.01 2.01 1.96 1.92

160 х 160 а£ 1.62 1.90 2.01 1.91 1.57 1.70

аа 1.57 1.69 2.03 2.03 1.85 1.84

160 х 320 а£ 1.84 1.78 1.94 2.01 1.99 1.68

аа 1.75 1.27 2.02 2.01 1.86 1.87

Таблица2

Параметр Функция

ф п Т и V

а£ аа —0.21 —0.60 0.15 0.51 1.02 0.67 0.28 0.47 0.50 0.49

описанной в предыдущих разделах, но с некоторыми упрощениями. Так был произведен расчет при Ке = 1, Ма = 0, Ог = 0, Еит = 0, де = 0, Qc = 0, Qf = 0, вязкость переменная, геометрия области такая же, как в следующем разделе. Результаты приведены в табл. 1.

Из табл. 1 видно, что почти все порядки сходимости лежат в пределах 1.5... 2.3 и большинство близки к 2 (теоретическое значение). Исключение составляет порядок аа сходимости функции П при = 160, 4М = 320. Возможной причиной является то, что коэффициент а1г+_1/2т при п =1 близок к 2 (см. формулы (27), (23)), а при п = 0 равен нулю (см. формулу (28)). Другой возможной причиной является влияние ошибок округления.

Был также произведен расчет порядков сходимости "в себе" при 4Ж = 4М = 120, Ог = 0, Еит = 0, Пс = 0, Qf = 0, вязкость постоянная, остальные параметры приняты такими, как в разд. 2, геометрия области — как в следующем разделе, т. е. расчет конвекции Марангони при реальном (т. е. очень большом) числе Рейнольдса. Результаты приведены в табл. 2. Отличие большинства порядков от 1 (теоретическое значение) может быть объяснено такими причинами, как недостаточное число разбиений, влияние ошибок округления и разрывность коэффициента аП+1/2т при переходе п от 0 к 1.

6. Пример расчета

Зададим конформное отображение квадрата 0 < х < 1, 0 < у < 1 на область, занимаемую расплавом, по формулам

, . с вЬ х , . с вт(п — 0.5 + у)

г(х,у) = —---—-т-, г (х,у) = —-

сЬх — сов(п — 0.5 + у)' ' сЬ х — сов(п — 0.5 + у)

где с выбирается из условия г(1, 0) = Дс. Тогда

Н (х,у) = —---.

сЬ х — сов(п — 0.5 + у)

Сетку на отрезке 0 < х < 1 построим следующим образом. Положим

, _ 3.75Ет 3.75Ет

-х1/2 = Ытт(а п\ , -1/2

NH(0,0)' ' 2NH(1,0)'

Тогда в плоскости переменных r, z максимумы соответствующих шагов сетки будут примерно в N/3.75 раз меньше толщины погранслоя в окрестности оси симметрии и свободной границы. (Для свободной границы роль толщины погранслоя играет величина em/2, для остальных границ она примерно в два раза больше.) Отношение /1/2 примем

постоянным и большим 1 при 1 < n < N/ 2 — 1. Это же отношение при N/2 + 1 < n < N — 1 примем постоянным и меньшим 1. Значения постоянных подберем так, чтобы hxN/2-1/2 = hxw/2+1/2 (здесь искомой величиной является также ж^/2). Описанными условиями сетка определяется однозначно. Сетку на отрезке 0 < y < 1 построим аналогично, приняв

hyi/2 = hxi/2 N/M, hyM-1/2 = hxi/2N/M.

На рис. 1, а и б показаны расчетные сетки в плоскости Огг размерности 30 х 30 и 120 х 120 соответственно.

Функцию пондеромоторной силы / (х) зададим в виде

/(х) = 0.5(1 + вт(пх)), 0 < х < 1.

Приводимые ниже примеры расчета относятся к нестационарным течениям. Задача решалась следующим образом. На сетке 30 х 30 задавались аналитически некоторое начальное приближение функций Ф, П, Т, удовлетворяющее граничным условиям на линиях фазового перехода, и начальное приближение, тождественно равное нулю, для функции Ж Если решалась задача с вращением, то задавалась угловая скорость вращения нижнего Пс4(£) и верхнего П^(£) цилиндров по закону

Пс,М*) = (1 — е-0'044)Пс>/, £ е [0, 250).

При £ > 250 полагалось Пс>/4(£) = Пс>/. Шаг т при всех расчетах был равен 10-3, а число итераций для функции тока Ь = 10, что обеспечивало для сетки 120 х 120 норму невязки в Ь2 ~ 10-6 — 10-7 при норме невязки начального приближения ~ 10-4. Сначала производили от 300 до 750 тысяч итераций по времени (в зависимости от варианта) на сетке

Рис. 1.

Рис. 2.

30 х 30. Затем все функции интерполировались на сетку размерности 30 х 60 и результат брался за начальное приближение при расчете на этой сетке. Производили 15 тысяч итераций и делали интерполяцию на следующую сетку, попеременно удваивая число узлов то по у, то по х, пока не достигалась размерность сетки 120 х 120. На этой сетке производили 500 тысяч итераций, что соответствует 100 с размерного времени. Последний этап занимал более одних суток машинного времени на ПЭВМ Реп1шш-4 с частотой 3 ГГц.

На рис 2, а, б, в показаны соответственно изолинии функции тока Ф, температуры Т и модифицированной азимутальной компоненты скорости Ш при решении задачи без каких-либо упрощений (вариант 1). Далее показаны изолинии функции тока Ф при постоянной вязкости (г) (вариант 2), при отсутствии термокапиллярной и пондеромоторной сил (д)

Таблица3

Вариант Функция хтт утт тт Хтах утах тах

1 Ф 0.57 0.11 —0.050 0.86 0.91 0.090

1 ^ху 0 1.0 0.997 2.80

1 Т 34, 0 0.9992 0.59 37, 5

1, 2, 3 W 1.0 1.0 —2.7 1.0 0.0 0.655

2 Ф 0.61 0.11 —0.049 0.88 0.91 0.093

2 ^ху 0 1.0 0.997 2.74

3 Ф 0.83 0.14 —0.020 0.91 0.41 0.051

3 ^ху 0 1.0 0.37 0.216

4 Ф 0.81 0.84 —0.125 0.65 0.37 0.199

4 ^ху 0 1.0 0.31 1.160

(вариант 3) и при отсутствии вращения (е) (вариант 4). Параметр 0 был равен 0.8, а Y = 0.99. Число итераций на сетке 30 х 30 составляло 750 тысяч для вариантов 1 и 2, 400 тысяч для варианта 3 (так как через некоторое время итерации для этого варианта разошлись) и 300 тысяч для варианта 4 (так как было получено стационарное решение). Характеристики течений приведены в табл. 3 (здесь обозначено vxy = Vu2 + v2, xmin, ymin — координаты x, y точки минимума; xmax, ymax — координаты x, y точки максимума соответствующей функции).

Заметим, что на сетке 120 х 120 все полученные решения являются нестационарными.

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

[1] Веретенцев В.А. Построение разностной сетки в области с криволинейными границами с помощью конформного отображения // Актуальные вопросы прикл. математики. М.: Изд-во МГУ, 1989. С. 88-93.

[2] Пивоваров Ю.В. О построении ортогональной разностной сетки в криволинейном четырехугольнике // Вычисл. технологии. 2003. Т. 8, № 5. С. 94-101.

[3] Булеев Н.И. Пространственная модель турбулентного обмена. М.: Наука, 1989. 344 с.

[4] Воеводин А.Ф., Юшкова Т.В. Численный метод решения начально-краевых задач для уравнений Навье — Стокса в замкнутых областях на основе метода расщепления // Сиб. журн. вычисл. математики. 1999. Т. 2, № 4. С. 321-332.

[5] Овчарова А.С. Метод расчета стационарных течений вязкой жидкости со свободной границей в переменных вихрь — функция тока // ПМТФ. 1998. Т. 39, № 2. С. 59-68.

[6] Зверев В.Г. Об одном итерационном алгоритме решения разностных эллиптических уравнений // Вычисл. технологии. 1999. Т. 4, № 1. С. 55-65.

[7] Muhlbauer A., Muiznieks A. et al. Interface shape, heat transfer and fluid flow in the floating zone growth of large silicon crystals with the needle-eye technique //J. Crystal Growth. 1995. Vol. 151. P. 66-79.

[8] Регель А.Р., Глазов В.М. Физические свойства электронных расплавов. М.: Наука, 1980. 295 с.

[9] Пивоваров Ю.В. Одномерная тепловая задача о бестигельной зонной плавке в быстропе-ременном магнитном поле // Динамика сплошной среды: Сб. науч. тр. / РАН. Сиб. отд-ние. Ин-т гидродинамики. 1996. Вып. 111. С. 100-108.

[10] Пивоваров Ю.В. Параметрический анализ задачи о бестигельной зонной плавке в магнитном поле // Динамика сплошной среды: Сб. науч. тр. / РАН. Сиб. отд-ние. Ин-т гидродинамики. 2000. Вып. 116. С. 142-147.

[11] Кочин Н.Е., Кивель И.А. и др. Теоретическая гидромеханика. Ч. II. М.: Гос. изд-во физ.-мат. лит., 1963. 727 с.

[12] Лойцянский Л.Г. Механика жидкости и газа. М.: Наука, 1973. 736 с.

[13] Андреев В.К., Капцов О.В. и др. Применение теоретико-групповых методов в гидродинамике. Новосибирск: Наука, 1994. 318 с.

[14] ВоЕводин А.Ф., Остапенко В.В. и др. Проблемы вычислительной математики. Новосибирск: Изд-во Сиб. отд-ния РАН, 1995. 154 с.

[15] Патанкар С. Численные методы решения задач теплообмена и динамики жидкости. М.: Энергоатомиздат, 1984. 150 с.

[16] Пивоваров Ю.В. Условия монотонности факторизованной разностной схемы для эволюционного уравнения с двумя пространственными переменными // Вычисл. технологии. 2001. Т. 6, № 4. С. 81-91.

[17] Самарский А.А., Николаев Е.С. Методы решения сеточных уравнений. М.: Наука, 1978. 592 с.

Поступила в редакцию 12 декабря 2004 г.