CC BY
33
Вычислительные технологии
Том 11, № 1, 2006
ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ КОНВЕКЦИИ В ПЛАВАЮЩЕЙ ЗОНЕ
Ю. В. ПИВОВАРОВ Институт гидродинамики СО РАН, Новосибирск, Россия
Hydrodynamical part of the problem on melting in the floating zone embedded in the magnetic field is solved using previously constructed sequence of the orthogonal difference grids and the computational method developed earlier. A comparison with the results of other authors is presented.
Введение
Для получения монокристаллов кремния большого радиуса (более 5 см) используется метод бестигельной зонной плавки в магнитном поле, который состоит в следующем. Верхняя (заготовка) и нижняя (выращиваемый монокристалл) части цилиндрического вертикального образца медленно движутся вниз и вращаются в противоположных направлениях. Часть нижней границы заготовки покрыта жидкой пленкой, остальная часть граничит с плавающей зоной, находящейся между заготовкой и монокристаллом и поддерживаемой в жидком состоянии неподвижным источником высокочастотного электромагнитного поля — индуктором. Токи, наводимые индуктором, сосредоточены в тонком скин-слое, примыкающем к свободной границе расплава. Они приводят к выделению джоулева тепла и создают пондеромоторную силу, направленную ортогонально свободной границе, экспоненциально убывающую при удалении от нее и являющуюся одним из источников конвекции в расплаве.
В работе [1] осесимметричная задача о бестигельной зонной плавке в магнитном поле решена численно в полной постановке с упрощающим предположением о малости конвективных членов по сравнению с вязкими в уравнении импульса в скин-слое. Это предположение позволяет снести пондеромоторную силу из уравнения импульса в граничное условие для вихря и в результате не слишком мельчить сетку в окрестности свободной границы.
В настоящей работе производится расчет движения расплава кремния в плавающей зоне методом, описанным в работе [2], на последовательности ортогональных разностных сеток, построенных в [3] (форма области в [3] бралась близкой к приведенной в [1]) с более сосредоточенным тепловыделением и, как следствие, в несколько раз большей пон-деромоторной силой, чем в [1]. Оценки показывают, что при этом конвективные и вязкие члены в уравнении импульса в скин-слое имеют одинаковый порядок. Это подтверждают и численные расчеты. Поэтому пондеромоторная сила из уравнения импульса в граничное условие для вихря не сносится. Скорость расплава вне скин-слоя получается на порядок
© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2006.
больше, чем в [1], а решение в зависимости от формы функции пондеромоторной силы либо близко к стационарному, либо носит колебательный характер с периодом, на порядок меньшим, чем в [1]. Просчитаны также случай отсутствия пондеромоторной и термокапиллярной сил, рассмотренный в [1], для которого получено удовлетворительное согласие с результатами расчетов, произведенных в [1], и случай, когда действуют пондеромотор-ная, термокапиллярная силы или сила плавучести отдельно от других сил при отсутствии вращения. В последних трех случаях решение выходит на стационарный режим.
1. Размерные параметры задачи
Плотность р и поверхностное натяжение а расплава будем считать линейными функциями температуры Т (переменность плотности учитывается в приближении Буссинеска). Кроме того, р терпит скачок при фазовом переходе. Кинематическую вязкость представим функцией вида
V(Т) = аи + Ъи/(Т - ои).
Ниже приведены названия, обозначения и значения материальных констант [1, 4]. Символ "т" означает жидкую, "в" — твердую фазы кремния. Коэффициент аи вычисляется по формуле
ау Vm Ъ-и/ (Тт Су).
Введем систему координат тОг следующим образом: ось Ог направим вверх вдоль оси образца, а ось От проведем горизонтально вправо через нижнюю точку трехфазного контакта.
Кроме перечисленных нам понадобятся следующие размерные параметры: толщина скин-слоя [5] £т = с/^/2пт= 2.85 • 10-2 см; радиус монокристалла тс = 5.1 см; координата т точки контакта расплава с жидкой пленкой (верхней правой угловой точки плавающей зоны) Tf = 1.333 см; скорость протягивания монокристалла vc = -5.52 • 10-3 см/с; скорость протягивания заготовки Vf = vc(тc/тf )2 = -8.085 • 10-2 см/с (здесь не учитывается, что большая часть массы поступает в расплав из жидкой пленки, а величина Vf
Название параметра
Символ
Значение
Круговая частота тока Постоянная Больцмана Коэффициент серости Температура плавления Значения р при Т = Тт
Значение V при Т = Тт -йр/йТ при Т > Тт -йа/йТ
Коэффициенты уравнения для V
Ускорение свободного падения Удельная теплоемкость Коэффициент теплопроводности Скорость света
Коэффициент электропроводности
ио а\
/Ьт Т
Тт рт
Рв
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 /с
принята такой, какой она была бы, если бы заготовка имела радиус г/); угловая скорость вращения верхней ш/ = -2.16 рад/с и нижней шс = 5.24 • 10-1 рад/с границ фазового перехода [1]; характерный размер в плавающей зоне I = 1.5 см; характерная скорость движения расплава г^ = 10 см/с [6]; характерный перегрев расплава ДТ = 50 К [1]; мощность индуктора N = 9.442 • 1010 эрг/с [6].
Предположим, что в плавающей зоне мощность выделяющегося джоулева тепла составляет 0.6№0 (а остальная часть мощности выделяется в жидкой пленке и твердых частях образца). Тогда максимум напряженности магнитного поля на свободной границе Гт плавающей зоны вычисляется по формуле
1/2
Но = | 8 • 0.6Жо/ | /2Шо£ту г(в)/(в/вт)^в
где в — безразмерная длина дуги на Гт, отсчитываемая от нижней точки Гт; г (в) — безразмерная координата г соответствующей точки Гт; вт — значение в в верхней точке Гт; /(ж) — функция, описывающая форму поверхностной плотности источников джоулева тепла на Гт, заданная на единичном отрезке и имеющая единичный максимум (см. разд. 3). Эта формула следует из формул (3), (5) работы [6], описывающих функцию поверхностной плотности источников джоулева тепла и связь между ней и мощностью индуктора.
Для 1-го и 2-го вариантов задания пондеромоторной силы (см. разд. 3) получим значения
Я01 = 287.7г1/2/(с • см1/2), Н02 = 315.6г1/2/(с • см1/2).
2. Безразмерные критерии подобия
Выберем в качестве масштабов длины, скорости, времени, вязкости и температуры соответственно /, г^, /Д>1, ит и ДТ. Тогда в задаче появятся следующие безразмерные параметры:
^оо дртДт/3
Ие = — = 4688, Ог = ^—— = 9.703 • 10
^т рт^т
5
Н 2 I Н 2 I Еит1 = --= 685.1, Еит2 = --= 824.4,
Лл ш0Н011 0 ооа п„ Ш0Н021 ц 1П
= ё-Л7г = 9.226, ^е2 = ц-= 11.10,
8пртСртг1ДТ 8пртСртг1ДТ
Ма = ^Т! = 2.895 • 105, Ре = = 54.37,
рт Лт
Ус = — = -5.520 • 10-4, V/ = — = -8.085 • 10-3, г1 г1
П/ = ^ = -3.240 • 10-1, Пс = — = 7.860 • 10-2, 5 = — = 0.9091,
г1 г1 рт
\та1Тт * = 0.7377, Т0 = ^ = 34, Лт ДТ ' 0 ДТ '
А = — = 0.6173, в
АТи„
= 0.2985, С =
ДТ
33.22,
Яс = Т
3.4, Я, = ^ = 0. ' I
Ет = у = 1.900 ■ 10-2.
Здесь Ив — число Рейнольдса; Ог — число Грасгофа; Еит1, Еит2 — магнитные числа Эйлера; Qe2 — отношения характерных интенсивностей джоулева тепловыделения и
конвективного теплопереноса; Ма — число Марангони; Ре — число Пекле; Ус, V, — безразмерные скорости протягивания; Пс, — безразмерные скорости вращения; Б — отношение плотностей; Б1 — число Био; Т0 — безразмерная температура плавления; AV, BV, СV коэффициенты уравнения для безразмерной вязкости
V
А + В/(Т - С);
Яс — безразмерный радиус монокристалла; Я, — безразмерная координата г точки контакта расплава с жидкой пленкой; Ет — безразмерная толщина скин-слоя.
Кроме того, при расчете используется параметр проскальзывания А1 = Ы/р^т, где к — коэффициент проскальзывания. Значение А1 полагается равным 5000.
3. Задание функции пондеромоторной силы
Функции магнитного давления, пондеромоторной силы на Гт, поверхностной (проинтегрированной по нормали к Гт) и объемной плотности источников джоулева тепла на Гт отличаются только постоянными множителями и поэтому описываются одной функцией f (ж), заданной на единичном отрезке и имеющей единичный максимум (см. формулы (7), (8) работы [5]). При вычислении формы Гт в [3] для задания магнитного давления использовалась функция f (ж), определенная формулой
f (ж) =
1 — сов
сЬ(во(0, 75ж — 0, 5)) — сов
1/8
54
4 3,
где в0 = 0.1, = 9 град. Эта функция принимается в качестве 1-го варианта задания пондеромоторной силы на Гт. Ее график — кривая 1 на рис. 1.
0.8 0.6 0.4 0.2 0
0.2 0.4 0.6 0.8 1
Ь
С
V
V
2
Рис. 1.
Таблица 1.
к 0 1 2 3 4 5 6 7 8
I 100 4.0 5.3 9.1 22.9 35.0 61.0 84.9 97.0 100
X 100 0.0 1.7 6.3 13.7 18.0 23.0 29.4 36.0 44.1
к 9 10 11 12 13 14 15 16 17
I 100 98.7 92.7 85.0 68.0 35.3 21.2 9.3 6.4 3.7
X 100 52.5 60.9 67.2 76.1 85.3 87.8 92.4 96.2 100
На рис. 2 работы [7] для радиуса кристалла 3.8 см приведены полученные расчетом форма свободной границы плавающей зоны и распределение по ней плотности электрического тока ^. Отсекая верхнюю часть свободной границы (чтобы получить ее форму, близкой к приведенной в [1]) и учитывая, что / пропорциональна получим функцию / (х), приведенную в табл. 1. Ей соответствует ломаная 2 на рис. 1.
Чтобы получить непрерывную функцию, разложим эту таблично заданную функцию по полиномам Лежандра Хп(х), удовлетворяющим условию нормировки
J ХП(х)^х =1, п = 0,Ж 0
Для этого определим коэффициенты Ьп по формулам
17
(хк+
2
Ьп = ^ (Хк+12 Хк) (/к) + /к+А^к+О), п = 0, Жл,
к=0
где к — номер колонки в табл. 1. Заметим, что нельзя брать число большим из-за малого числа разбиений отрезка 0 < х < 1. Расчет показывает, что при = 7 кривая
/ (х) = ь„хх„(х)
п=0
хорошо аппроксимирует таблично заданную функцию внутри отрезка 0 < Х < 1, но вблизи его границ появляются участки немонотонности /(х). Чтобы их ликвидировать, преобразуем переменную х в правой части предыдущего равенства:
/ (х) = ^ Ь„Х„(0.97(х + 0.05)/1.05).
п=0
Сделаем также растяжение вдоль оси ординат, обеспечивающее максимум /(х), равный единице. Так определенная функция / (х) принимается в качестве 2-го варианта задания пондеромоторной силы на Гт. Ее график — кривая 3 на рис. 1.
В полуполосе 0 < в < вт, 0 < п< то функции пондеромоторной силы и объемной плотности источников джоулева тепла с точностью до постоянных множителей определяются формулой
/п(в,п) = /(в/вт) ехр(—2п/Ет),
где п — безразмерная внутренняя нормаль к Гт (см. формулу (5), следующую за ней формулу и формулу (18) работы [2]).
1
4. Результаты расчета конвекции
Конвекция расплава в плавающей зоне обусловлена такими факторами, как:
1) пондеромоторная сила (результат действия магнитного поля);
2) термокапиллярный эффект (зависимость коэффициента поверхностного натяжения от температуры);
3) сила плавучести, обусловленная зависимостью плотности расплава от температуры;
4) вращение образца;
5) протягивание образца.
Расчет производился методом, описанным в [2], на последовательности ортогональных разностных сеток, построенных в [3]. При расчете задавалось шесть различных вариантов исходных данных. В пяти вариантах, за исключением 4-го, функции пондеромоторной силы и объемной плотности источников джоулева тепла определялись в соответствии с 2-м вариантом задания функции f (ж) (см. предыдущий разд.), а в 4-м варианте — в соответствии с 1-м вариантом задания функции f (ж). Протягивание образца и переменность вязкости учитывались во всех шести вариантах расчета. 1-й, 2-й и 3-й варианты расчета соответствуют действию первого, второго, третьего факторов, обусловливающих конвекцию, отдельно от других факторов. В 4-м и 5-м вариантах учитывалось действие всех факторов, а в 6-м варианте — действие 3-, 4- и 5-го факторов.
Опишем порядок расчетов. Пусть жО'у — правая система координат в плоскости независимых пространственных переменных задачи, такая, что функции г (ж, у), г (ж, у) конформно отображают прямоугольник П: 0 < ж < 1, 0 < у < У = 1.08732, на двумерную область, занятую расплавом (в работе [3] константа У обозначалась через /), так что отрезки у = 0, у = У, ж = 0, ж = 1 переходят соответственно во фронт кристаллизации Гс, фронт плавления Г,, ось симметрии Г0 и свободную границу Гт.
Начальные приближения для функции тока Ф и модифицированного вихря П брались при расчете удовлетворяющими граничным условиям при у = 0, у = У и линейными по у. (Для Ф это обеспечивало и выполнение граничных условий при ж = 0 и ж =1, а для П — граничного условия при ж = 0.) Начальное приближение для температуры Т бралось удовлетворяющим граничным условиям при у = 0 и у = У и параболическим по у с максимумом, равным Т0 + 1 (это обеспечивает перегрев расплава 50 К, граничное условие при ж = 0 для Т при этом также выполняется). Если решалась задача с вращением, то начальное приближение для модифицированной азимутальной компоненты скорости Ш бралось равным нулю, а угловые скорости вращения нижнего и верхнего цилиндров задавались зависящими от времени £ по закону
Псл^) = (1 — е-4)Пс;/, £ е [0,10).
При £ > 10 полагалось Пс>,4(£) = Пс,,. Шаг по времени т во всех расчетах был равен 10-3, что соответствует 1.5 ■ 10-4 с размерного времени, малость т обусловлена малостью минимального шага сетки в переменных ж, у : ж^ — ж^-1 = 4.275 ■ 10-3/Ж, где N — число разбиений по оси О'ж (см. [3, с. 96]), весовой параметр 7 = 0.99, а коэффициент нижней релаксации при решении задачи для Ф 0 = 0.8.
Число итераций Ь при решении задачи для Ф на каждом шаге по времени определялось следующим образом. Задавалось некоторое малое число . Производилось Ьт;п итераций.
Если
—
\
Ф^Н2аЫу/ | ^ Ф2я^5 +
дП П
/0 Н
<е0ф,
где Фп — функция невязки, Н — (дг/дж)2 + (дг/ду)2 ция правой части при решении задачи для Ф, дП —
— коэффициент Ламэ, /ф — функ-граница П, в — длина дуги дП,
то итерации прекращались. В противном случае они продолжались до тех пор, пока еф не станет меньше еф или число итераций не станет равным Ьтах. Число итераций К по времени при решении всей задачи определялось аналогично числу итераций при решении
задачи для Ф с заменой Ьт
Ьтах на Kmin,
Кт
еф — на еП и е0 — на
еп —
\
(дП/д£)2Н2^Ыу/ // П2Н2^Ыу.
П
Для всех вариантов было Ьт^ — 2, Ьтах — 20, еф — еП — 10-8. Сначала задача для вариантов 1-3 решалась на сетке 30 х 30 с Кт^ — 104, Ктах — 5 • 105, затем — на сетках 30 х 60, 60 х 60 и 60 х 120 каждый раз с начальным приближением, построенным на предыдущей сетке, при Ктщ — Ктах — 1.5 • 104. И окончательный расчет делался на сетке 120 х 120 при КтЫ — 104, Ктах — 5 • 105.
Для вариантов 1-3 число итераций по времени составило на сетке 30 х 30 соответственно 76 499, 199 966 и 344 737, а на сетке 120 х 120 — 67 185, 358 165 и 500 000. Значение параметра еп в 3-м варианте на сетке 120 х 120 составило 1.2 • 10-8. На последнем шаге по времени число итераций при решении задачи для Ф у вариантов 1-3 на сетках 30 х 30 и 120 х 120 составляло Ь — 2. Значение еф начального приближения Ф на последнем шаге по времени на сетке 30 х 30 для вариантов 1-3 составило 1.3 • 10-11, 4.4 • 10-12 и 3.2 • 10-12, а на сетке 120 х 120 — 1.2 • 10-11, 1.5 • 10-11 и 2.9 • 10-13 соответственно. После проведения двух итераций значение еф на сетке 30 х 30 для вариантов 1-3 составило соответственно 9.9 • 10-16, 5.6 • 10-15 и 7.7 • 10-15, а на сетке 120 х 120 — 3.0 • 10-13, 1.0 • 10-12 и 1.5 • 10-13.
Попытки получить численное решение по такой же схеме для четвертого варианта оказались неудачными: на сетках 30 х 30 и 30 х 60 процесс расчета заканчивался аварийным остановом. Поэтому для четвертого — шестого вариантов задача решалась следующим образом: производилось 6 • 104 итераций по времени на сетке 60 х 60, затем 104 итераций на сетке 60 х 120 и окончательно 5 • 105 итераций на сетке 120 х 120. На последнем шаге по времени значение еф начального приближения при решении задачи для Ф для вариантов 4-6 на сетке 120 х 120 составило соответственно 2.5 • 10-4, 1.9 • 10-9 и 2.4• 10-7. Число итераций Ь при решении задачи для Ф составило для четвертого — шестого вариантов 20, 2 и 10 соответственно, значение еф по окончании решения задачи для Ф — 1.1 • 10-7, 3.7 • 10-10 и 9.1 • 10-9 соответственно, а величина еп по окончании расчетов — 5.9, 2.3 • 10-4 и 2.0 • 10-1 соответственно.
На рис. 2 показаны изолинии функции тока (а — е) для первого — шестого вариантов соответственно, а на рис. 3 — изотермы для этих же вариантов. На рис. 4, а — в показаны изолинии модифицированной азимутальной компоненты скорости Ш соответственно для четвертого — шестого вариантов.
В табл. 2 приведены характеристики течений для шести вариантов расчета (здесь угг — модуль проекции вектора скорости расплава на плоскость гОг; гтщ, и гтах, гтах —
О 0.5 1 1.5 2 2.5 3 Г
0 0.5 1 1.5 2 2.5 3 г
0 0.5 1 1.5 2 2.5 3 г
Рис. 2.
0 0.5 1 1.5 2 2.5 3 г
0 0.5 1 1.5
0 0.5 1 1.5 2 2.5 3 г
Рис. 3.
0 0.5 1 1.5
0 0.5 1 1.5 2 2.5 3 г
0 0.5 1 1.5 2 2.5 3
Рис. 4.
Таблица 2.
Функция Вариант гтт -2тт тт гтах ^тах тах
ф 1 2.94 0.043 -0.233 1.82 -0.175 1.154
2 2.77 0.030 -0.177 1.60 -0.255 0.095
3 2.65 -0.034 -0.075 0.770 -0.076 0.034
4 2.90 0.033 -0.321 1.82 -0.175 1.120
5 2.95 0.015 -0.193 1.74 -0.189 1.003
6 2.70 -0.018 -0.070 1.13 -0.076 0.037
УГ2 1 — — 0 1.56 0.396 3.781
2 — — 0 3.40 0.002 3.021
3 — — 0 0.199 -0.263 0.244
4 — — 0 0 0.058 9.479
5 — — 0 1.56 0.396 3.615
6 — — 0 0.343 -0.134 0.220
Т 1 0.885 0.891 -0.239 1.03 0.558 39.51
2 — — 0 1.29 0.446 85.61
3 — — 0 2.57 0.395 122.0
4 0.001 -1.138 -9.38 0.880 0.764 40.33
5 0.887 0.898 -0.072 1.20 0.476 42.85
6 — — 0 2.02 0.384 126.8
W 4, 5, 6 0.888 0.902 -0.256 3.40 0 0.909
координаты г, г точек минимума и максимума соответственно; единицей измерения скорости служит г^ — 10 см/с, а температуры — градус Кельвина; началом отсчета температуры является точка плавления).
Из табл. 2 видно, что в первом, четвертом и пятом вариантах расчета минимальная температура расплава получилась меньше нуля. Но сама зона отрицательных температур очень мала. В первом и пятом вариантах она примыкает к границам Гт, Г/, а в четвертом варианте — к границам Г0, Гс. Анализ выходного массива Т показывает, что в плоскости гОг эта зона помещается в квадрат со стороной 3 • 10-2 см, 5.4 • 10-3 см и 4.5 • 10-3 см соответственно для первого, четвертого и пятого вариантов расчета. Существование такой зоны может быть объяснено тем, что форма границы плавающей зоны взята близкой к приведенной в [1], а условия задачи в настоящей работе существенно отличны от принятых в [1] (главное различие — более концентрированное тепловыделение и как следствие — большая пондеромоторная сила).
На рис. 5, а показано распределение касательной скорости на свободной границе (в — безразмерная длина дуги, отсчитываемая от нижней точки свободной границы, положительным для скорости является направление в сторону возрастания в) и на различных расстояниях от нее для пятого варианта расчетов. Кривая 1 — распределение скорости на самой свободной границе, 2 — на расстоянии Ет, 3 — на расстоянии 2Ет и 4 — на расстоянии 3Ет от нее. На рис. 5, б показано распределение скорости вдоль внутренней нормали п, за единицу длины которой принята величина Ет, при в — 1.5.
На рис. 6, а приведено распределение температуры Т на свободной границе и на различных расстояниях от нее также для пятого варианта расчетов. Кривая 1 — распределение температуры на самой свободной границе, 2 — на расстоянии 3Ет, 3 — на расстоянии 6Ет и 4 — на расстоянии 9Ет от нее, рис. 6, б — аналог рис. 5, б с заменой функции касательной скорости на функцию температуры.
а б
Рис. 5.
аб
Рис. 6.
Наконец, на рис. 7 приведены зависимости минимума — кривые 1, 3, 5 и максимума — кривые 2, 4, 6 функции тока от размерного времени соответственно для четвертого, пятого и шестого вариантов расчета.
Проанализируем полученные результаты. Движение, осуществляющееся в первом варианте расчетов, можно назвать электроконвекцией, так как оно вызвано только силой электрической природы (влияние протягивания образца ввиду малой его скорости несущественно). Максимум скорости в этом варианте достигается вблизи свободной границы и составляет около 38 см/с. Характерная скорость внутри области течения (т.е. за границами пограничных слоев для скорости) составляет 10... 20 см/с. Этого оказывается достаточно для формирования температурных пограничных слоев (сгущения изотерм вблизи границ расплава — см. рис. 3, а). Скорость на свободной границе (аналог кривой 1 на рис. 5, а) не имеет отрицательного всплеска вблизи точки в — 0 и положительного — вблизи точки в — 3 (эти всплески у кривой 1 на рис. 5, а обусловлены термокапиллярным эффектом).
Движение, осуществляющееся во втором варианте расчетов, можно назвать термокапиллярной конвекцией. Максимум скорости в этом движении составляет 30 см/с и достигается на свободной границе. Но характерная скорость внутри области на порядок меньше.
1.6 ^^
1.2
2 4
0.8
0.4
0
Г
б 5 3
1
-0.4
О
20
40
60
Рис. 7.
Поэтому температурных пограничных слоев вблизи границ расплава не образуется.
Движение, осуществляющееся в третьем варианте расчетов, принято называть термогравитационной конвекцией. Максимум скорости достигается вблизи оси симметрии и составляет 2.4 см/с. Так как движение очень медленное, тепловых пограничных слоев в этом варианте также не образуется.
Таким образом, предположение, сделанное в [6], о том, что характерная скорость внутри области течения составляет примерно половину от максимальной скорости в пограничном слое вблизи свободной границы, оказалось выполненным для случая электроконвекции и не выполненным для случая термокапиллярной конвекции. Выводы работы [6] о том, что электроконвекция доминирует и что характерная скорость расплава составляет ~ 10 см/с, оказались верными. Толщину возникающих динамического е.и и теплового £т пограничных слоев можно теоретически оценить по формулам
Сравнение теоретических , £т и фактических е'Г1), е'т (см. рис. 5, б и 6, б) толщин пограничных слоев показывает, что £-ю ^ , £т ~ 2еТ, т.е. согласие между ними вполне приемлемое.
Движение, осуществляющееся в четвертом и пятом вариантах расчета, можно назвать комбинированной конвекцией, так как оно обусловлено действием всех факторов в совокупности. Но между ними есть качественное различие, вызванное разным видом функции пондеромоторной силы. В четвертом варианте на оси симметрии возникает струя с максимумом скорости 95 см/с. Радиус участка, где скорость превышает 30 см/с, составляет 1 мм, а его высота — около 1 см (как показали дополнительные расчеты, при отсутствии вращения и прочих равных условиях этого участка в окрестности оси симметрии нет). Течение
в четвертом варианте расчетов носит колебательный характер, о чем говорят высокое значение величины еп = 5.9 и вид кривых 1, 2 на рис. 7, иллюстрирующих зависимость минимума и максимума функции тока от времени для этого варианта. Период колебаний обеих величин составляет 0.2025 с, временной интервал вывода их на график — 0.0075 с. (Под периодом здесь понимается расстояние между соседними максимумами, движение при этом может и не быть периодическим.) Размах колебаний составляет 0.0255 для минимума и 0.0244 для максимума функции тока, а среднее по периоду колебаний значение равно -0.3092 для минимума и 1.1315 для максимума функции тока. В пятом варианте расчета максимум скорости на оси симметрии составляет 2.6 см/с, т.е. движение вблизи оси симметрии очень медленное. Глобальный максимум скорости достигается вблизи свободной границы и составляет 36 см/с. Невысокое значение величины еп = 2.3 • 10-4 по окончании расчетов в пятом варианте и вид кривых 3, 4 на рис. 7, иллюстрирующих зависимость минимума и максимума функции тока от времени для этого варианта, говорят о том, что решение в этом случае вышло на режим, близкий к стационарному. Начиная с момента t = 6.03 с до момента t = 75 с значение минимума функции тока монотонно возрастает от -0.20275 до -0.19273, а значение максимума функции тока при 4.92 c < t < 75 c монотонно убывает от 1.03515 до 1.00265.
Движение, осуществляющееся в шестом варианте расчетов, назовем термогравитационной конвекцией с вращением. Максимум скорости в меридиональной плоскости достигается для этого варианта вблизи оси симметрии и составляет 2.2 см/с. Ввиду медленности движения тепловых пограничных слоев вблизи границ области, занятой расплавом, не образуется. Течение является нестационарным и носит колебательный характер, о чем говорят довольно высокое значение величины еп = 0.20 и вид кривых 5, 6 на рис. 7. Минимум функции тока испытывает колебания, не видные на рисунке, ввиду малости их размаха, равного 0.002. Размах колебаний максимума функции тока составляет 0.018. Период колебаний в обоих случаях равен 6 с. Средние по периоду колебаний значения составляют -0.070695 для минимума и 0.040965 для максимума функции тока.
Вид рис. 2, г — е и 4, а — в показывает, что модифицированная азимутальная компонента скорости W претерпевает изменения, необходимые для удовлетворения граничных условий, только в зонах малых скоростей. В зонах больших скоростей она близка к константе (для каждой зоны своя константа). Это происходит потому, что функция W = const
2W fdvdrdvdr\
с точностью до малых членов---— — —--+ — — (ввиду гармоничности функции
rH2Re \ dx dx dy dy J
r(x,y) и формул (3) работы [2]) удовлетворяет уравнению (14) работы [2]. Это явление аналогично явлению постоянства вихря ш в плоской задаче в зоне больших скоростей, детально разобранному в работе [8], и явлению постоянства температуры внутри области течения в первом, четвертом и пятом вариантах, менее выраженному ввиду относительно небольшого числа Пекле. Такое же явление должно иметь место и для модифицированного вихря П в первом и втором вариантах расчета, в которых отсутствуют сила плавучести и вращение.
Сравним результаты расчетов с полученными в [1]. Из способов обезразмеривания, примененных в настоящей работе и работе [1], следует, что функцию тока ф работы [1] нужно преобразовать в функцию Ф для сравнения с функцией Ф настоящей работы следующим образом:
Ф = - ^Ф = -0.00077ф,
V\i2
где = 3.4 • 10-3 см2/с — коэффициент кинематической вязкости, используемый в [1].
Первый и второй варианты в [1] соответствуют комбинированной конвекции с разбиением области течения на 7760 и 3006 треугольных элементов соответственно. Третий вариант соответствует термогравитационной конвекции с вращением. Приводимые ниже данные являются универсальными для всех трех вариантов, рассчитанных в [1]. Средние значения минимума и максимума функции ф
Ф min = -0.058, Ф max = 0.037,
максимальная скорость (ее можно оценить из рис. 7 [1]) vmax = 0.25, т.е. 2.5 см/с. Перегрев расплава в среднем 42.5 K. Период колебаний минимума и максимума функции ф — примерно 4 и 7 с соответственно. Размах колебаний — 0.0024 для минимума и 0.012... 0.16 для максимума функции ф.
Таким образом, результаты расчета, сделанного в настоящей работе и в [1] для комбинированной конвекции, отличаются принципиально (за исключением максимума температуры, для которого имеет место случайное совпадение), а для случая термогравитационной конвекции с вращением — довольно близки (за исключением максимума температуры). Последний факт можно объяснить следующим образом. Хотя перегрев расплава в настоящей работе больше полученного в [1] в три раза, но используемый коэффициент рт ~ в 2.3 раза меньше, так что эффективное число Грасгофа больше только на 30%. Эффективное число Пекле в обеих работах примерно одинаково, коэффициент теплопроводности (а следовательно, и эффективное число Био) в настоящей работе меньше на 22 %, вязкость — на 5-30 %, используется другая функция источников джоулева тепла, не учитывается тепловое излучение индуктора и используются несколько другие граничные условия для вихря на твердых стенках и свободной границе. Таким образом, различия в постановке задачи хотя и имеют место, но не очень велики. Это и обеспечивает близость результатов расчета термогравитационной конвекции с вращением в обеих работах, в том числе и близость изолиний функции тока (см. рис. 2, е настоящей работы и рис. 7, d работы [1]). Разницу в максимальной температуре (перегреве) расплава для рассматриваемого варианта можно объяснить только разной долей выделяющегося в плавающей зоне джоулева тепла (60 % в настоящей работе и в несколько раз меньшей долей — в работе [1]).
Автор выражает благодарность В. В. Кузнецову за полезное обсуждение результатов работы.
Список литературы
[1] 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.
[2] Пивоваров Ю.В. Расчет движения жидкости с переменной вязкостью в области с криволинейной границей // Вычисл. технологии. 2005. Т. 10, № 3. С. 87-107.
[3] Пивоваров Ю.В. О построении ортогональной разностной сетки в криволинейном четырехугольнике // Вычисл. технологии. 2003. Т. 8, № 5. С. 94-101.
[4] Регель А.Р., Глазов В.М. Физические свойства электронных расплавов. М.: Наука. 1980. 295 с.
[5] Пивоваров Ю.В. Одномерная тепловая задача о бестигельной зонной плавке в быстропе-ременном магнитном поле // Динамика сплошной среды: Сб. науч. тр. / РАН. Сиб. отд-ние. Ин-т гидродинамики. 1996. Вып. 111. С. 100-108.
[6] Пивоваров Ю.В. Параметрический анализ задачи о бестигельной зонной плавке в магнитном поле // Динамика сплошной среды: Сб. науч. тр. / РАН. Сиб. отд-ние. Ин-т гидродинамики. 2000. Вып. 116. С. 142-147.
[7] Lie K.H., Walker J.S. ет al. Melt motion in the float zone process with an axial magnetic field // J. Crystal Growth. 1991. Vol. 109. P. 167-173.
[8] Betchelor G. K. On steady laminar flow with closed streamlines at large Reynolds number // J. Fluid Mech. 1956. Vol. 1, N 2. P. 177-190.
Поступила в редакцию 26 августа 2005 г., в переработанном виде — 24 ноября 2005 г.
