авторефераты диссертаций БЕСПЛАТНАЯ БИБЛИОТЕКА РОССИИ

КОНФЕРЕНЦИИ, КНИГИ, ПОСОБИЯ, НАУЧНЫЕ ИЗДАНИЯ

<< ГЛАВНАЯ
АГРОИНЖЕНЕРИЯ
АСТРОНОМИЯ
БЕЗОПАСНОСТЬ
БИОЛОГИЯ
ЗЕМЛЯ
ИНФОРМАТИКА
ИСКУССТВОВЕДЕНИЕ
ИСТОРИЯ
КУЛЬТУРОЛОГИЯ
МАШИНОСТРОЕНИЕ
МЕДИЦИНА
МЕТАЛЛУРГИЯ
МЕХАНИКА
ПЕДАГОГИКА
ПОЛИТИКА
ПРИБОРОСТРОЕНИЕ
ПРОДОВОЛЬСТВИЕ
ПСИХОЛОГИЯ
РАДИОТЕХНИКА
СЕЛЬСКОЕ ХОЗЯЙСТВО
СОЦИОЛОГИЯ
СТРОИТЕЛЬСТВО
ТЕХНИЧЕСКИЕ НАУКИ
ТРАНСПОРТ
ФАРМАЦЕВТИКА
ФИЗИКА
ФИЗИОЛОГИЯ
ФИЛОЛОГИЯ
ФИЛОСОФИЯ
ХИМИЯ
ЭКОНОМИКА
ЭЛЕКТРОТЕХНИКА
ЭНЕРГЕТИКА
ЮРИСПРУДЕНЦИЯ
ЯЗЫКОЗНАНИЕ
РАЗНОЕ
КОНТАКТЫ


Pages:   || 2 |
-- [ Страница 1 ] --

РОССИЙСКАЯ АКАДЕМИЯ НАУК

СИБИРСКОЕ ОТДЕЛЕНИЕ

ИНСТИТУТ ГИДРОДИНАМИКИ им. М.А. ЛАВРЕНТЬЕВА

На правах рукописи

Пивоваров Юрий Владимирович

МОДЕЛИРОВАНИЕ КОНВЕКЦИИ РАСПЛАВА

ПОЛУПРОВОДНИКОВОГО МАТЕРИАЛА ПРИ

ЗОННОЙ ПЛАВКЕ

05.13.18 математическое моделирование, численные методы и

комплексы программ

Диссертация на соискание ученой степени кандидата физико-математических наук

Научный руководитель:

Доктор физ.-мат. наук, Воеводин А. Ф.

Новосибирск – 2006 СОДЕРЖАНИЕ ВВЕДЕНИЕ 4 ГЛАВА 1. ПОСТРОЕНИЕ ОРТОГОНАЛЬНОЙ РАЗНОСТНОЙ СЕТКИ 1.1. Построение одномерной сетки, сгущающейся на краях отрезка 1.2. Вычисление формы свободной границы 1.3. Построение неортогональной двумерной сетки 1.4. Ортогонализация разностной сетки ГЛАВА 2. МОДЕЛИРОВАНИЕ КОНВЕКЦИИ ПРИ БЕСТИГЕЛЬНОЙ ЗОННОЙ ПЛАВКЕ В МАГНИТНОМ ПОЛЕ 2.1. Размерные параметры задачи 2.2. Безразмерные критерии подобия 2.3. Электромагнитная часть задачи 2.4. Гидродинамическая часть задачи 2.5. Численный алгоритм 2.6. Тестовые расчеты 2.7. Результаты расчета конвекции ГЛАВА 3. УСЛОВИЯ МОНОТОННОСТИ ФАКТОРИЗОВАННОЙ РАЗНОСТНОЙ СХЕМЫ ДЛЯ ПАРАБОЛИЧЕСКОГО УРАВНЕНИЯ 3.1. Описание разностной схемы 3.2. Достаточные условия монотонности в общем случае 3.3. Необходимые и достаточные условия монотонности для задачи с малым числом разбиений 3.4. Пример расчета ГЛАВА 4. ЗАДАЧА О ВРАЩАЮЩЕМСЯ КОНТЕЙНЕРЕ 4.1. Постановка задачи 4.2. Асимптотика вихря и функции тока в окрестности точек контакта 4.3. Аналитические решения 4.4. Метод и результаты расчета течения ЗАКЛЮЧЕНИЕ СПИСОК ЛИТЕРАТУРЫ ВВЕДЕНИЕ Зонная плавка кристаллографический метод рафинирования материалов, состоящий в перемещении узкой расплавленной зоны вдоль длинного твер дого стержня из рафинируемого материала. В настоящей диссертации рас сматривается два варианта зонной плавки.

1 вариант: бестигельная зонная плавка (БЗП) в магнитном поле (МП), применяемая, в частности, для выращивания монокристаллов кремния боль шого радиуса (более 5 см). Верхняя (заготовка) и нижняя (выращивае мый монокристалл) части вертикального цилиндрического образца медлен но движутся вниз и вращаются в противоположных направлениях с разными угловыми скоростями. Часть нижней границы заготовки покрыта жидкой пленкой, остальная часть граничит с плавающей зоной, находящейся между заготовкой и монокристаллом, поддерживаемой в жидком состоянии непо движным источником высокочастотного электромагнитного поля индук тором, и удерживаемой между твердыми частями образца силами поверх ностного натяжения и магнитного давления. Токи, наводимые индуктором, сосредоточены в тонком скин-слое, примыкающем к свободной границе рас плава. Они приводят к выделению джоулева тепла и создают пондеромотор ную силу, направленную ортогонально свободной границе, экспоненциально убывающую при удалении от нее, и являющуюся одним из источников кон векции в расплаве.

Полученные методом БЗП в МП монокристаллы кремния большого ра диуса используются в основном в двух направлениях:

1) в силовой электронике создание тиристоров, силовых транзисторов, используемых в мощных силовых преобразователях;

2) при изготовлении высокоэффективных солнечных батарей.

2 вариант: зонная очистка поликристаллического полупроводникового материала (в частности, германия). Цилиндрический полый контейнер на полнен полупроводниковым материалом и расположен под небольшим углом к горизонту. Он вращается и одновременно совершает медленное поступа тельное движение вдоль своей оси, при этом часть его нагревается до высо кой температуры. Материал полупроводника плавится. Образуются фронты плавления и кристаллизации. Расплав не целиком заполняет контейнер: име ется свободная поверхность. Предполагается также, что расплав отделяет от стенки контейнера тонкий слой мелкодисперсной смазки.

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

Первое упоминание о методе зонной плавки относится к 1927 г., когда он был использован П. И. Капицей для получения однородных монокристал лов висмута и Маккеханом для очистки железа [9]. Широкую известность он получил благодаря работам В. Дж. Пфанна (США), который применил его в 1952 г. для получения германия высокой степени чистоты в контей нере (не вращающемся). Метод БЗП был впервые предложен П. Кек и М.

Голей (США) в 1953 г. для очистки полупроводникового кремния. В первых установках радиус очищаемого образца составлял 0.2...1 см, так что произ водительность установок была весьма низкой. Нынче этот радиус составляет уже 5...10 см.

Важным является вопрос о распределении температуры в образце при БЗП, так как от него зависят размер и форма расплавленной зоны. Д. До нальд [13] в 1961 г. предложил одномерную модель зонной плавки, обла дающую аналитическим решением. Н. Кобаяши [49] численно исследовал осесимметричную задачу в цилиндрической области. В [2] численно рассчи таны поле температуры и форма свободной границы для радиуса образца 0.75 см в условиях невесомости. Во всех названных выше работах гидроди намическое течение расплава не учитывалось. Однако в [5] показано, что оно существенно влияет на величину и форму жидкой зоны. Причем в за висимости от температуры на торцах образца оно может приводить как к уменьшению, так и к увеличению ее размеров.

В [52] численно решена осесимметричная стационарная задача о БЗП в земных условиях в полной постановке (найдены форма границы плавающей зоны, поля скоростей и температуры) без учета вращения, а в [51] с уче том вращения образца для его радиуса 0.2 см. Однако в этих двух работах зависящая от температуры вязкость введена в уравнения для вихря и ази мутальной компоненты скорости неправильно (ср. формулы (1), (2) [51] и формулы (2.4.4), (2.4.5), (2.4.14) при x = r, y = z, H = 1 настоящей диссертации).

Во всех упомянутых выше работах нагрев образца производился путем его облучения.

В [43] изучен характер особенности переменного во времени электромаг нитного поля вблизи вершины проводящего клина, моделирующего угловую точку плавающей зоны.

В [41, 42] теоретически рассмотрена задача об определении границы фронта плавления и свободной границы в окрестности точки трехфазно го контакта применительно к задаче о БЗП в МП. Действие индуктора при этом заменяется поверхностным тепловыделением и магнитным давлением.

В [40] решена аналогичная задача с учетом термокапиллярных сил.

В [25, 27] рассмотрена и численно решена модельная стационарная осе симметричная задача о БЗП в МП при отсутствии вращения. Однако маг нитное поле индуктора здесь учитывалось только как источник внешнего давления на свободную поверхность и не учитывалось как источник элек троконвекции в расплаве.

В [56] численно решена осесимметричная нестационарная задача о БЗП в МП в земных условиях с учетом вращения для радиуса монокристалла 5.15 см в полной постановке, но с некоторыми упрощающими допущениями.

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

Действие этой формулы сносится на свободную границу плавающей зоны и она используется в качестве граничного условия для вихря. Таким образом, течение внутри скин-слоя исключается из рассмотрения, в результате чего отпадает необходимость сильно мельчить сетку вблизи свободной границы расплава. В [55] рассчитывается распределение примеси в монокристалле кремния, получаемом методом БЗП в МП. При этом поле скоростей в пла вающей зоне рассчитывается методом, описанным в [56], и осредняется по времени, а затем решается стационарная задача о распределении примеси в расплаве и растущем монокристалле. В [58] решается задача о БЗП в МП и рассчитывается распределение примеси в монокристалле уже в трехмерной постановке. Наконец, в работах [57] и [50] решались задачи, аналогичные решенным в [55] и [58] соответственно, но при наличии дополнительного низкочастотного индуктора.

Однако анализ показывает, что указанное допущение, делаемое в послед них пяти из названных работ, на практике не выполняется. Поэтому акту альной является разработка модели процесса БЗП в МП с учетом пондеро моторной силы не в граничном условии, а в самом уравнении для вихря. При этом сетка в окрестности свободной границы расплава должна быть значи тельно более мелкой, чем сетка, используемая в [56]. Из-за этого должен быть существенно уменьшен и временной шаг, что влечет увеличение числа итераций по времени на порядок. Но с учетом увеличения производитель ности персональных электронно-вычислительных машин со времени выхода работы [56], время счета этой задачи в современных условиях сравнимо со временем счета задачи, решенной в [56] более 10 лет назад, и составляет 1– суток.

При разработке модели процесса БЗП встает вопрос о методе решения задачи.

В [55, 56, 57] использовался метод конечных элементов с разбиением об ласти на треугольные ячейки и со сгущением расчетной сетки вблизи всех границ плавающей зоны.

В [25, 27] использовался конечно-разностный метод сквозного счета на неортогональной сетке, причем преобразование координат состоит в простом растяжении радиальной переменной, так что исходной областью независи мых переменных становится прямоугольник (см. [26]). В [51, 52] использован конечно-разностный метод на неортогональной сетке общего вида, так что области, занятые заготовкой, плавающей зоной и выращиваемым монокри сталлом, отображаются на прямоугольники, причем на границах фазового перехода все узлы сетки являются общими для каждой пары граничащих об ластей. В обоих случаях сетки в исходной прямоугольной области являются равномерными, но если в [25, 27] сетка в физической области также близка к равномерной, то в [51, 52] она сгущается вблизи всех границ плавающей зоны, так как там возможно образование пограничных слоев для скорости.

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

В настоящей диссертации для решения гидродинамической части осе симметричной нестационарной задачи о БЗП в МП впервые применяется конечно-разностный метод с использованием консервативной монотонной разностной схемы [7, с. 280]. Обычно в таких задачах применяют либо схему с центральными разностями для конвективных членов, являющуюся кон сервативной, но обладающую свойством монотонности только при малых скоростях течения [26] (при больших скоростях применение этой схемы при водит к пилообразному решению и может вызвать расходимость итераций по времени), либо схему с разностями против потока, являющуюся монотон ной при любых скоростях движения, но не обладающую свойством консер вативности [56]. Схема из [7] совмещает достоинства обеих схем, содержит первую как предельный случай при скорости, стремящейся к нулю, и близ ка ко второй при скорости, стремящейся к бесконечности. Кроме того, в настоящей диссертации впервые при решении задачи о БЗП применяется ортогональная сетка, получающаяся с помощью конформного отображения прямоугольника на область, занятую расплавом [8]. Сетка в прямоугольнике строится сгущающейся на его границах таким образом, чтобы в произволь ном нормальном сечении к границе физической области внутрь погранс лоя, образующегося около этой границы, попало-бы не менее 6 точек при максимальном числе разбиений области. Преимущества ортогональной раз ностной сетки очевидны: в левых частях уравнений отсутствуют смешанные производные. Ясно, что сгустить ее вблизи линии раздела вихрей нельзя, но использование консервативной схемы позволяет надеяться на то, что, хотя погранслой на линии раздела вихрей будет размываться схемной вязкостью, значения искомых функций вне его будут стремиться к правильным пре делам (то есть к решению дифференциальной задачи) [28]. При решении разностных уравнений используется метод раздельного решения уравнений для вихря и функции тока, описанный в [10, 26]. Впервые в этой задаче правильно учтена в уравнениях для вихря и азимутальной компоненты ско рости зависящая от температуры вязкость и, самое главное, в уравнении для вихря учтена объемная пондеромоторная сила, сосредоточенная внутри скин-слоя.

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

В [3, 14, 39, 44, 46, 47] изучались в основном режимы движения, когда жидкость полностью покрывает стенки цилиндра, а внутри имеется газо вая полость. Так, в [14, 47] аналитическими методами исследовалось пове дение слоя жидкости в быстро вращающемся цилиндре. При быстром вра щении центробежные силы доминируют над капиллярными и гравитаци онными (действие капиллярных сил в [14, 47] не учитывалось). Свободная поверхность жидкости близка к цилиндрической. Экспериментальное изуче ние такого режима течения проводилось в [46]. В экспериментальной работе [44] при малых скоростях вращения обнаружены трехмерные стационарные течения. Они возникают при увеличении скорости от нуля еще до того, как жидкость полностью покроет боковую поверхность цилиндра. В [3] получе ны некоторые необходимые условия существования и достаточные условия несуществования плоских течений в двусвязной области в медленно враща ющемся цилиндре. В [39] решалась задача о ползущем движении жидкости во вращающемся цилиндре методом граничных элементов в отсутствие ка пиллярных сил также для случая, когда жидкость полностью покрывает стенки цилиндра.

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

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

Диссертация состоит из введения, четырех глав, заключения и списка литературы. Главы разделены на разделы, а разделы четвертой главы на пункты. Нумерация формул выполнена раздельно для каждого раздела и остается сквозной внутри него. Нумерация рисунков и таблиц единая для всей диссертации. Список литературы составлен по алфавиту русскому и английскому отдельно. Работы автора приведены в конце списка литерату ры.

В первой главе строится последовательность ортогональных разностных сеток применительно к задаче о БЗП в МП методом, описанным в [8]. Форма области при этом берется близкой к рассчитанной в [56]. (Она удовлетвори тельно согласуется с данными экспериментальных измерений.) Начальное приближение неортогональная сетка строится методом трансфинитной интерполяции [20, 21]. Во второй главе осуществляется постановка задачи, описывается численный алгоритм и производится численный расчет гидро динамического течения в расплаве при БЗП в МП на последовательности ортогональных разностных сеток, построенных в главе 1. В третьей главе выводятся условия монотонности разностной схемы, используемой при ре шении задачи о БЗП в МП. В четвертой главе аналитически и численно исследуется задача о движении расплава во вращающемся контейнере при малых и умеренных скоростях его вращения.

Настоящая диссертация выполнена в Институте гидродинамики им. М.

А. Лаврентьева СО РАН, Новосибирск, под руководством А. Ф. Воеводи на. Основные результаты диссертации опубликованы в 5 главе монографии [61], в статьях в сборниках [62, 64] и в статьях в журналах [63, 65, 66, 67] и докладывались соавторами работ на Всесоюзном семинаре по математи ческому моделированию процессов кристаллизации (Рига, 1989 г.), а также автором на десятой Зимней Школе по механике сплошных сред (Пермь, 1995 г.), на третьей Сибирской школе-семинаре "Математические проблемы механики сплошных сред" (Новосибирск, 1999 г.) и Всероссийской конферен ции "Теория и приложения задач со свободными границами" (Бийск, г.). Кроме того, они были доложены на семинаре отдела прикладной гид родинамики Института гидродинамики им. М. А. Лаврентьева СО РАН и на объединенном семинаре "Информационно–вычислительные технологии" под руководством академика Ю. И. Шокина и профессора В. М. Ковени в Институте вычислительных технологий СО РАН. Материалы диссертации использовались в двух отчетах ИГиЛ СО РАН.

В заключение этого введения автор выражает благодарность В. В. Пухна чеву за постановку задач, научному руководителю А. Ф. Воеводину, а также В. В. Кузнецову, О. М. Лаврентьевой, А. С. Овчаровой и В. И. Яковлеву за помощь при выполнении работы.

ГЛАВА ПОСТРОЕНИЕ ОРТОГОНАЛЬНОЙ РАЗНОСТНОЙ СЕТКИ В этой главе решается задача о построении ортогональной разностной сет ки в области, имеющей форму криволинейного четырехугольника, приме нительно к задаче о БЗП в МП. Вычисляется форма свободной границы плавающей зоны при некоторой заданной аналитически функции магнитно го давления. Задается форма фронтов плавления и кристаллизации, близкая к приведенной в [56]. В полученной области строится неортогональная сетка по формулам трансфинитной интерполяции [20, 21]. Ортогонализация сетки производится по алгоритму, описанному в [8]. Основные результаты главы изложены в работе автора [63].

1.1 Построение одномерной сетки, сгущающейся на краях отрезка Пусть N четное число;

h1min, h2min заданные вещественные положи тельные числа, меньшие 1/N. Определим на отрезке 0 x 1 сетку xn, n = 0, N, следующим образом:

n h1min · (z1 1)/(z1 1), при 0 n N/2, xn = (1.1.1) 1 h2min · (z2 n 1)/(z2 1), при N/2 n N, N где z1, z2 неизвестные пока вещественные числа, бльшие единицы. По о ложим z1 = 1 + 2k1 /N, z2 = 1 + 2k2 /N (1.1.2) и обозначим через = xN/2 (1.1.3) точку стыковок сеток при n N/2 и n N/2. Из условий стыковки следует:

k1 = ((1 + 2k1 /N )N/2 1)N/(2), (1.1.4) k2 = ((1 + 2k2 /N )N/2 1)N/(2(1 )). (1.1.5) Эти уравнения при заданном служат для определения k1, k2. Они реша ются методом деления отрезка пополам. После этого z1 и z2 восстанавли ваются по формулам (1.1.2). Для определения наложим дополнительное условие:

xN/2+1 xN/2 = xN/2 xN/21 = hmax. (1.1.6) Из определения hmax, z1 и z2 следует, что N/ 1 1 1 1 1/z hmax + +... + N/21 = hmax =, 1 z1 1 1/z z N/ 1 1 1 1 1/z hmax + +... + N/21 = hmax = 1.

1 z2 1 1/z z Исключая отсюда путем деления одного на другое неизвестную hmax и выражая, получим:

=. (1.1.7) N/ (11/z1 )(11/z2 ) 1+ N/ (11/z2 )(11/z1 ) Так как z1 и z2 здесь зависят от в силу (1.1.2), (1.1.4), (1.1.5), то (1.1.7) представляет собой нелинейное уравнение для определения. Оно легко решается итерационно: сначала полагаем = 0.5 и находим значения z1, z2. Затем подставляем их в (1.1.7) и вычисляем новое значение и т.

д.

В численных расчетах при каждом решении уравнений (1.1.4), (1.1.5) ис пользуется 100 итераций, при решении уравнения (1.1.7) такое-же коли чество итераций. Итого значения функций правых частей (1.1.4), (1.1.5) вы числяются 10000 раз.

После решения уравнений (1.1.7), (1.1.4), (1.1.5) с учетом (1.1.2) сетка строится по формулам (1.1.1). В качестве исходных данных для ее построе ния нужно задать N, h1min, h2min.

1.2. Вычисление формы свободной границы Пусть s длина дуги, отсчитываемая от нижней точки свободной границы m плавающей зоны, Hs | амплитуда колебаний во времени s-й компо ненты напряженности МП на = m l, l граница жидкой пленки, q поверхностная плотность источников джоулева тепла;

fn функция магнитного давления, N0 = 9.442·1010 эрг/с мощность индуктора, вычис ленная для радиуса образца 5 см, 0 = 1.76 · 107 рад/с круговая частота тока в индукторе, m = 2.85 · 102 см толщина скин-слоя в жидкой фазе.

Имеют место соотношения (см. раздел 2.3):

q (s) = (0 m /(16))(Hs | )2, (1.2.1) fn (s) = (1/(16))(Hs | )2, (1.2.2) N0 2r(s)q (s)ds. (1.2.3) Последнее равенство имеет место при условии, что джоулевым тепловыде лением в твердой фазе кремния можно пренебречь.

Введем функцию f (x), заданную на единичном отрезке и имеющую единичный максимум, описывающую форму функций магнитного давления и джоулева тепловыделения:

fn (s) = (H0 /(16))f (s/sm ), (1.2.4) где H0 максимум Hs |m, s безразмерная длина дуги, sm безраз мерная длина кривой m. В качестве масштаба длины выберем l = 1.5 см.

Предположим, что мощность токов, выделяющаяся на m, равна 0.6N0.

Тогда на основании (1.2.1)–(1.2.4) для H0 получим следующее выражение:

sm H0 = (8 · 0, 6N0 /(l2 0 m r(s)f (s/sm )ds))1/2. (1.2.5) Система уравнений для функций r(s), z(s), параметрически задающих границу m, в безразмерных переменных имеет вид [56, 59]:

r (s) = sin (s), z (s) = cos (s), cos (s) (s) = + Bo(z(s) + F ef (s/sm ) c0 ), (1.2.6) r(s) где m gl2 H Bo =, Fe =, (1.2.7) m 16m gl m = 2.53г/см3 плотность жидкого кремния, g = 980см/с2 ускоре ние свободного падения, m = 737 дин/см коэффициент поверхностного натяжения расплава.

К (1.2.6) нужно присоединить условия Коши:

r(0) = Rc, z(0) = Zc, (0) = c, (1.2.8) где c = 11 град угол между касательной к свободной границе при s = и направлением z, при котором процесс БЗП протекает устойчиво (выра щиваемый кристалл имеет постоянный радиус [60]).

Для определения констант c0, sm необходимо поставить 2 дополнитель ных условия, например, r(sm ) = Rf ;

z(sm ) = Zf. (1.2.9) Функция f (x) должна, по видимому, резко возрастать от значения, близ кого к нулю, в окрестности точки x = 0, и затем плавно меняться до точки x = 1. Зададим эту функцию в виде 1/ 1 cos f (x) = ·54 · 4 3, ch(0 (0, 75x 0, 5)) cos где 0 = 0.1;

0 = 9 град. Ее график показан на рис. 1.

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

Сначала задача (1.2.6) – (1.2.9) решалась при H0 = 316.9 Гс (оценочное f 0. 0. 0. 0. x 0 0.2 0.4 0.6 0.8 Рис. значение H0 ). Затем по формуле (1.2.5) было определено H0 = 287.7 Гс.

После этого был произведен окончательный расчет формы свободной гра ницы. Задача (1.2.6) – (1.2.9) решалась методом Эйлера с помощью визу альной пристрелки по параметру c0. (Для каждого значения c0 стро ился график свободной границы, после чего значение c0 подправлялось по методу деления отрезка пополам.) Исходные данные были следующие:

Rc = 3.4, Zc = 0, c = 11 град, Rf = 0.8884, Zf = 0.9020. Число раз биений по s было принято равным 256. Результаты расчетов следующие:

F e = 0.4428;

Bo = 7.569;

sm = 3;

c0 = 0.801898. График свободной грани цы приведен в следующем разделе (рис. 3).

1.3. Построение неортогональной двумерной сетки Пусть в плоскости rOz дан криволинейный четырехугольник ABCD (рис. 2), причем границы AB и CD однозначно проектируются на ось Or и заданы уравнениями z D C A B O r Рис. 0 AB : z = f11 (r), r11 r r12, 0 CD : z = f21 (r), r21 r r22, а границы BC и DA однозначно проектируются на ось Oz и заданы уравнениями 0 BC : r = f12 (z), z11 z z12, 0 DA : r = f22 (z), z21 z z22.

Пусть в плоскости xO y дан единичный квадрат K : 0 x 1, 0 y 1, в котором задана неравномерная сетка 0 = x0 x1... xN 1 xN = 1, 0 = y0 y1... yM 1 yM = 1. Построим взаимно-однозначное отображение K на ABCD. Будем говорить, что сетка n, введенная на отрезке [a, b], подобна сетке xn, если n = a(1 xn ) + bxn, n = 0, N.

Аналогично будем говорить, что сетка m, введенная на отрезке [c, d], подобна сетке ym, если m = c(1 ym ) + dym, m = 0, M.

0 0 0 0 1 Введем на отрезках [r11, r12 ] и [r21, r22 ] сетки r1n и r2n, подобные 0 0 0 0 1 сетке xn, а на отрезках [z11, z12 ] и [z21, z22 ] сетки z1m и z2m, подобные сетке ym, и положим согласно методу трансфинитной интерполяции (см.

формулы (2.4), (2.5) работы [21]) 1 rnm = f12 (z1m )(1 xn ) + f22 (z2m )xn, 1 znm = f11 (r1n )(1 ym ) + f21 (r2n )ym, n = 0, N, m = 0, M.

Тогда если 1 min(f22 (z2m ) f12 (z1m )) m и 1 min(f21 (r2n ) f11 (r1n )) 0, n то множество пар (rnm, znm ) устанавливает взаимно однозначное соот ветствие между точками x = xn, y = y m квадрата K и точками r = rnm, z = znm криволинейного четырехугольника ABCD, причем точки, лежащие на сторонах K, переходят в точки, лежащие на соответ ствующих частях границы ABCD.

Построим криволинейный четырехугольник, форма которого является типичной для плавающей зоны (см. [56]). Зададим фронт кристаллизации (линию AB) уравнением z = Zc 1.14(1 r2 /Rc ), 0 r Rc, фронт плавления (линию CD) уравнением z = Zf 0.27(1 + cos(180r/Rf )), 0 r Rf и ось симметрии (линию BC) уравнением r = 0, Zc 1.14 z Zf 0.54, где Rc = 3.4, Zc = 0, Rf = 0.8884, Zf = 0.9020 безразмерные геометри ческие параметры (в качестве единицы измерения служит величина l = 1. см). Форма свободной границы (линии AD) была вычислена в предыдущем разделе. Так как AD не проектируется однозначно на ось Oz, повернем построенный криволинейный четырехугольник по часовой стрелке на град, построим в нем сетку, как это описывалось выше, и сделаем обратный поворот. При построении сетки в ABCD нужно задать сетку в квадрате K. Она строилась сгущающейся на границах K таким образом, чтобы максимальные шаги уже ортогонализованной сетки в ABCD (см. следу ющий раздел) в направлении нормали к соответствующим границам были примерно в N/15 (M/15) раз меньше толщины погранслоя в окрестности соответствующей границы (для границы AD эта толщина составляет 0. см, для остальных границ она примерно в 2 раза больше). Это приводит к следующим соотношениям:

h12min = xN xN 1 = 4.275 · 103 /N, h11min = x1 x0 = 16 · h12min, h21min = y1 y0 = 8.55 · 103 /M, h22min = yM yM 1 = 16 · h21min.

Способ построения одномерных сеток на отрезках 0 x 1, 0 y изложен в разделе 1.1. На рис. 3 показана неортогональная сетка, построен ная по описанному выше алгоритму, при N = M = 30.

1.4 Ортогонализация разностной сетки Для построения ортогональной сетки в области D0, тождественной криво линейному четырехугольнику ABCD, используем алгоритм, основанный на минимизации функционала 2 2 2 1 x y 1 x y (x, y, Y ) = +Y + +Y drdz, Y r r Y z z D описанный в [8]. Здесь x(r, z), y(r, z) функции, отображающие D0 на K, Y числовая переменная.

Рис. Уравнениями Эйлера для функций x(r, z), y(r, z), выражающими необ ходимое условие существования экстремума функционала, являются урав нения Лапласа, которые после обращения роли зависимых и независимых переменных и деления на (((r/y)2 + (z/y)2 )((r/x)2 + (z/x)2 ))1/ принимают вид:

2r 2r 2r 2z 2z 2z a 2 2b + c 2 = 0, a 2 2b + c 2 = 0, (1.4.1) x xy y x xy y где (r/y)2 + (z/y)2 a=, c=, (r/x)2 + (z/x)2 a (r/x)(r/y) + (z/x)(z/y) b=.

((r/y)2 + (z/y)2 )((r/x)2 + (z/x)2 ) Условие трансверсальности на границе области D0 после обращения роли зависимых и независимых переменных принимает вид r r z z + = 0, (1.4.2) x y x y (r(0, y), z(0, y)) BC, (r(1, y), z(1, y)) DA, (r(x, 0), z(x, 0)) AB, (r(x, 1), z(x, 1)) CD. (1.4.3) Оно выражает ортогональность линий x = const и y = const на границе области D0.

На решении задачи a Y ;

b 0, так что уравнения (1.4.1) после замены x = x;

y = Y y переходят в уравнения Лапласа. Функции r(, y ) = x r(, y /Y ) и z (, y ) = z(, y /Y ) осуществляют конформное отображение x x x прямоугольника : 0 x 1;

0 y Y на D0.

Опишем итерационный процесс построения ортогональной разностной сетки в D0. Сначала при фиксированном расположении узлов сетки на границе определяются положения внутренних узлов путем приближенного решения разностных аналогов уравнений (1.4.1). Это приближенное реше ние рассчитывается попеременно - треугольным методом при фиксирован ных функциях a, b, c [34].

Далее при фиксированных внутренних узлах ищутся новые положения граничных узлов. Для этого аппроксимируем (1.4.2). Например, для грани цы BC можем написать k1 k k k k (x1 x0 )c2 r0m+1 r0m r0m r0m r1m r0m 0m + x1 x0 (ym+1 ym1 ) ym+1 ym ym ym k1 k1 k (x1 x0 )c r0m+1 r0m1 z1m z0m 0m + + (ym+1 ym1 ) x1 x0 (ym+1 ym1 ) k1 k k1 k k k z0m+1 z0m z0m+1 z0m z0m z0m = 0, (1.4.4) ym+1 ym ym ym1 (ym+1 ym1 ) где k номер итерации. Это равенство имеет второй порядок аппрокси мации на решении при условии ym+1 2ym + ym1 = O((ym+1 ym )2 ).

k k Чтобы найти r0m, z0m, нужно определить пересечение прямой (1.4.4) с гра ницей BC. Каждая из четырех границ D0 задается уравнениями вида r = ri (s), z = zi (s), где ri (s), zi (s) параболические сплайны переменной s [36]. Так что для нахождения очередного приближения положений точек на границе нужно последовательно решать некоторые квадратные или ли Рис. нейные (как в случае границы BC) уравнения относительно s, пока не будет найден корень, принадлежащий данному промежутку.

Итак, ортогональная разностная сетка строится в результате многократ ного чередования описанных итерационных процессов. Заметим, что для ее построения не требуется однозначной проектируемости участков границы на соответствующие оси.

На рис. 4 показана ортогональная сетка в D0 размерности 120 120.

Она была построена следующим образом. Сначала была ортогонализована сетка размерности 30 30, построенная в предыдущем разделе. Путем ин терполяции эта сетка была преобразована к размерности 30 60 и вновь ортогонализована. Процесс удвоения числа узлов сетки с последующей ее ортогонализацией продолжался, пока не была построена сетка максималь ной размерности. Число итераций при решении уравнения (1.4.2) с учетом (1.4.3) было равно 10. Обозначим через P число итераций при решении уравнений (1.4.1) попеременно–треугольным методом;

Q число чередо ваний итерационных процессов построения внутренних и граничных узлов сетки;

(r/x)(z/y) (r/y)(z/x) sin = ((r/y)2 + (z/y)2 )((r/x)2 + (z/x)2 ) синус угла между координатными линиями;

(Y r/x z/y) · 100% (r/y + Y z/x) · 100% S1 =, S2 = (Y r/x)2 + (r/y)2 (Y r/x)2 + (r/y) относительные погрешности в выполнении условий Коши - Римана, где Y определяется 1 раз после завершения итерационного процесса построения сетки размерности 120 120. Введем также интегральные характеристики:

S = sin dxdy;

1 = S1 dxdy;

K K 2 b2 dxdy · 100%.

2 = S2 dxdy;

b = K K В табл. 1 приведены значения параметров P, Q, N, M, S, 1, 2, b при построении ортогональной сетки. Здесь первая строка соответствует сет ке, построенной в предыдущем разделе. Заметим, что большие значения па раметра Q на начальном этапе построения сетки размерности 3030 необ ходимы ввиду очень малого шага сетки на границе K по нормали к ней.

Суммарное время расчета сетки размерности 3030 на ПЭВМ Pentium – 3 с частотой 1.2 Ггц составляет 8.5 минут. (При необходимости его можно суще ственно уменьшить, если на начальном этапе временно ввести равномерную сетку в K.) При P 128 время счета можно приближенно определить по формуле tc = 5 · 106 P QN M с.

Из табл. 1 видно, что при удвоении числа узлов (в каждом из направлений O x, O y ) сетки размерности 30 30 она значительно приближается к ортогональной (параметры S и b ), а при удвоении числа узлов сетки размерности 6060 этого уже не происходит. Вероятной причиной является влияние ошибок округления.

Таблица P Q N M S 1 2 b 0 0 30 30 0.88213 50.3 36.6 – 4 4000 30 30 0.98945 10.2 11.7 14. 8 2000 30 30 0.99915 3.3 4.0 4. 6 1000 30 30 0.99968 2.4 2.6 2. 32 500 30 30 0.99977 2.2 2.4 2. 64 250 30 30 0.99980 2.2 2.3 2. 128 125 30 30 0.99981 2.2 2.3 1. 128 125 30 60 0.99989 1.8 2.0 1. 128 60 60 60 0.999971 1.4 1.0 0. 128 60 60 120 0.999978 1.2 1.0 0. 128 60 120 120 0.999976 1.1 0.9 0. 256 60 120 120 0.999977 1.1 0.9 0. 512 1 120 120 0.999938 1.1 1.0 512 1 60 120 0.999946 1.1 1.1 512 1 60 60 0.999941 1.1 1.1 512 1 30 60 0.999905 1.4 1.8 512 1 30 30 0.99988 1.7 2.1 Заключительным этапом построения ортогональной сетки в D0 является определение Y по формуле Y= adxdy, K решение уравнений (1) при a Y, b 0, c 1/Y с граничными условиями 1-го рода (см. последние 5 строк таблицы) и растяжение переменной y в Y раз. Значение Y = 1.08732 определено по сетке размерности 120 и затем приписано решению на более редких сетках.

Сравнение результатов расчетов при различных N иM показало, что решение имеет 1.5-й порядок сходимости, хотя все уравнения и гранич ные условия аппроксимированы со 2-м порядком. Это объясняется нали чием степенной особенности в тупом угле A области D0. Анализ пове дения в окрестности точки A локальных характеристик, соответствующих 12-й строке табл. 1, показывает следующее. Параметр sin терпит мини мум, равный 0.9984, а параметр b = b · 100% максимум модуля, равный 5.7%, при n = 119, m = 1 (x = 0.99996, y = 7.1 · 105 ). (Все харак теристики определяются только внутри D0.) При n 115, m 2 (x 0.99977, y 1.5 · 104 ) имеем sin 0.9997, |b | 2.5%, а вне этого прямоугольника sin 0.9997, |b | 2.5%. То есть параметры sin, b, отвечающие за ортогональность сетки, несколько ухудшаются лишь в малой окрестности точки A. Ухудшение параметров S1 и S2, отвечающих за равенство коэффициентов Ламэ в направлениях x, y, более значительно и сосредоточено на бльших площадях: max |S1 | = 22.6% достигается при о n = 118, m = 1 (x = 0.99992, y = 7.1 · 105 );

|S1 | в основном больше 6% при n 76, m 28 (x 0.9441, y 0.01177);

max |S2 | = 9.0% дости гается при n = 119, m = 32 (x = 0.99996, y = 0.018342);

|S2 | частично больше 6% при n 92, m 40 (x 0.9921, y 0.04371).

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

ГЛАВА МОДЕЛИРОВАНИЕ КОНВЕКЦИИ ПРИ БЕСТИГЕЛЬНОЙ ЗОННОЙ ПЛАВКЕ В МАГНИТНОМ ПОЛЕ В этой главе производится расчет осесимметричной нестационарной кон векции жидкого кремния в плавающей зоне при заданной форме ее грани цы на построенной ранее (см. главу 1) последовательности ортогональных разностных сеток. Задача ставится в переменных вихрь – функция тока в соответствии с приближением Буссинеска. Дифференциальные уравнения для функции тока, вихря, азимутальной компоненты скорости и темпера туры записываются в дивергентной форме. Дифференциальные операторы конвективного и диффузионного переноса аппроксимируются монотонны ми консервативными разностными операторами, имеющими второй порядок аппроксимации по пространственным переменным при малых числах Рей нольдса и Пекле и первый при больших числах Рейнольдса и Пекле (здесь используются результаты работы [7, с. 280–288]). При решении разностных уравнений используется метод, позволяющий точным образом разделить за дачи вычисления вихря и функции тока, описанный в работах [10, 26]. При решении на каждом временнм шаге уравнения для функции тока исполь о зуется эффективный итерационный алгоритм, не требующий информации о свойствах и границах спектра разностного оператора, обобщающий ал горитм, описанный в работе [15]. Расчет производится для двух вариантов задания пондеромоторной силы на свободной границе и для них получают ся принципиально разные режимы конвекции. В первом варианте максимум и минимум функции тока колеблются с периодом 0.2 с и движение суще ственно нестационарно. Во втором варианте максимум и минимум функции тока изменяются со временем монотонно и движение постепенно выходит на режим, близкий к стационарному. Основные результаты главы изложены в работах автора [62, 64, 65, 67].

2.1. Размерные параметры задачи Плотность и поверхностное натяжение расплава будем считать линей ными функциями температуры T (переменность плотности учитывается в приближении Буссинеска). Кроме того, терпит скачек при фазовом переходе. Кинематическую вязкость представим функцией вида (T ) = a + b /(T c ).

В табл. 2 приведены названия, обозначения и значения материальных кон стант [32, 56]. Символ ”m” означает жидкую, ”s” твердую фазу кремния.

Коэффициент a вычисляется по формуле a = m b /(Tm c ).

Пусть Orz цилиндрическая система координат, где r полярный радиус, полярный угол. В сечении = const ось Or проведем го ризонтально вправо через нижнюю точку трехфазного контакта, а ось Oz вертикально вверх вдоль оси образца. Кроме перечисленных нам понадо бятся следующие размерные параметры: толщина скин-слоя 2m 0 = 2.85 · 102 см, m = c/ радиус монокристалла rc = 5.1 см, координата r точки контакта рас плава с жидкой пленкой (верхней правой угловой точки плавающей зоны) rf = 1.333 см, скорость протягивания монокристалла vc = 5.52 · 103 см/с, скорость протягивания заготовки vf = vc (rc /rf )2 = 8.085·102 см/с (здесь не учитывается, что бльшая часть массы поступает в расплав из жидкой о пленки, а величина vf принята такой, какой она была-бы, если-бы заготовка имела радиус rf ), угловая скорость вращения верхней f = 2.16 рад/с и нижней c = 5.24 · 101 рад/с границ фазового перехода [56], характерный размер в плавающей зоне (равный половине ее высоты) l = 1.5 см, характер ная скорость движения расплава v1 = 10 см/с (см. раздел 2.7), характерный Таблица Название параметра Символ Значение 1.76 · 107 рад/с Круговая частота тока 5.67 · 105 эрг/(с· см2 K4 ) Постоянная Больцмана Коэффициенты серости fbm 0. fbs 0. Температура плавления Tm 1700 K 2.53 г/см Значения при T = Tm m 2.30 г/см s 3.2 · 103 см2 /с Значение при T = Tm m 1.52 · 104 г/(см3 К) d/dT при T Tm T d/dT T 0.1дин/(см·К) 4.776 · 102 K·см2 /c Коэффициенты ур-я для b c 1661 K 980 см/с Ускорение своб. падения g 7.45 · 106 см2 /(К·c2 ) Удельная теплоемкость cpm 5.2 · 106 г·см/(К· с3 ) К-нты теплопроводности m 5.2 · 106 г·см/(К· с3 ) s 3 · 1010 см/с Скорость света c 1016 1/c К-нт электропроводности m перегрев расплава T = 50 K [56], мощность индуктора N0 = 9.442 · эрг/с (см. раздел 2.3).

Максимальная на m амплитуда колебаний во времени касательной к m компоненты вектора напряженности магнитного поля при сделанных в разделе 2.3 предположениях определяется по формуле sm H0 = (8 · 0.6N0 /(l2 0 m r(s)f (s/sm )ds))1/2, где s безразмерная длина дуги на m, отсчитываемая от нижней точки m, r(s) безразмерная координата r соответствующей точки m, sm значение s в верхней точке m, f (x) функция, описывающая форму поверхностной плотности источников джоулева тепла на m, заданная на единичном отрезке и имеющая единичный максимум. Для 1-го и 2-го вари антов задания f (x) (см. раздел 2.3) получим значения H01 = 287.7 Гс, H02 = 315.6 Гс.

2.2. Безразмерные критерии подобия Выберем в качестве масштабов длины, скорости, времени, вязкости и тем пературы соответственно l, v1, l/v1, m и T. Тогда в задаче появятся следующие безразмерные параметры:

gT T l v1 l = 9.703 · 105, Re = = 4688, Gr = m m m 2 H01 l H02 l Eum1 = = 685.1, Eum2 = 2 = 824.4, 8m m v1 8m m v 2 0 H01 l 0 H02 l Qe1 = = 9.226, Qe2 = = 11.10, 8m cpm v1 T 8m cpm v1 T T T l v1 lm cpm = 2.895 · 105, P e = Ma = = 54.37, m m m vc vf = 5.520 · 104, Vf = = 8.085 · 103, Vc = v1 v f l c l s = 3.240 · 101, c = = 7.860 · 102, S = f = = 0.9091, v1 v1 m fbm 1 Tm l Tm Bi = = 0.7377, T0 = = 34, m T T a b c A = = 0.6173, B = = 0.2985, C = = 33.22, m T m T rc rf m = 1.900 · 102.

Rc = = 3.4, Rf = = 0.8884, Em = l l l Здесь Re число Рейнольдса, Gr число Грасгофа, Eum1, Eum2 маг нитные числа Эйлера, Qe1, Qe2 отношения характерных интенсивностей джоулева тепловыделения и конвективного теплопереноса, M a число Ма рангони, P e число Пекле, Vc, Vf безразмерные скорости протягивания, c, f безразмерные скорости вращения, S отношение плотностей, Bi число Био, T0 безразмерная температура плавления, A, B, C коэффициенты уравнения для безразмерной вязкости = A + B /(T C ), Rc безразмерный радиус монокристалла, Rf безразмерная координата r точки контакта расплава с жидкой пленкой, Em безразмерная толщина скин-слоя.

2.3. Электромагнитная часть задачи Если 0, 0, V0, L0 характерные значения проводимости, промежутка времени, скорости и длины, а c скорость света, то при выполнении условий V 1 V 1, 1, 2 1 (2.3.1) 0 0 0 L0 c в уравнениях Максвелла можно пренебречь токами смещения (1/(4)) E/t, а в законе Ома конвективными токами e v (см.

[19, с. 45]), где e плотность заряда, v скорость среды, E вектор напряженности электрического поля, t время. (Впрочем, конвективные токи в нашей задаче не возникают, т. к. e = 0 в силу ее осевой симметрии.) В задаче о БЗП в МП 0 = m = 1016 c1 (проводимость жидкого кремния), 0 = 108 c (период колебаний тока), V0 = 40 см/с (максимальная скорость расплава в зоне действия электромагнитного поля), L0 = 3 см. Так что условия (2.3.1) выполнены с большой точностью. Кроме того, в силу малости магнитного числа Рейнольдса 40 V0 L = 1.7 · Rem = c не будем учитывать в законе Ома движение среды, т. е. отбросим член (m /c)v H, где H вектор напряженности магнитного поля. Тогда урав нения Максвелла примут вид (см.[19, с. 15] или [24, с. 21]):

rotH = j, (2.3.2) c 1 H rotE =, (2.3.3) c t div E = 4e, (2.3.4) div H = 0. (2.3.5) А закон Ома вид (см. [19, с. 45] или [24, с. 21]):

j = E, (2.3.6) где проводимость среды кусочно постоянная функция. Из (2.3.2), (2.3.6) следует:

rotH = E. (2.3.7) c Из (2.3.5) H = rotA, (2.3.8) где A некоторая вектор–функция, называемая векторным потенциалом.

Из (2.3.3), (2.3.8) 1 A rot(E + ) = 0. (2.3.9) c t Откуда 1 A + E = grad, (2.3.10) c t где функция, называемая скалярным потенциалом [16, с. 207].

Из (2.3.7), (2.3.8), (2.3.10) c A + rotrotA = grad. (2.3.11) t Положим div A = 0.

Так как rotrotA = A + graddiv A (2.3.13) (см. [16, с. 207]), то c A A = grad. (2.3.14) t В осесимметричном случае ненулевыми являются компоненты A, H r, H z, j, E. Предположим, что все функции совершают гар монические колебания во времени с круговой частотой 0. Введем в j рассмотрение комплексные амплитуды A, Hr, Hz,, E,, полагая A (r, z, t) = Re[A (r, z)ei0 t ], Hr (r, z, t) = Re[Hr (r, z)ei0 t ], Hz (r, z, t) = Re[Hz (r, z)ei0 t ], j (r, z, t) = Re[ (r, z)ei0 t ], j E (r, z, t) = Re[E (r, z)ei0 t ], (r, z,, t) = Re[(r, z, )ei0 t ]. (2.3.15) (Здесь у функции допускается линейная зависимость от, что не нарушает осевую симметрию задачи.) Тогда для A получим уравнение 2 A + c 1 (rA ) c i0 A + =. (2.3.16) z 4 r r r r На линиях разрыва следует поставить условия [A ] = 0, [ A /n] = 0, (2.3.17) где квадратные скобки означают скачек заключенной в них функции при переходе из одной среды в другую. Эти условия необходимы и достаточны для непрерывности вектора H при переходе через границу раздела двух сред.

На оси симметрии для A должно выполняться условие симметрии A = 0 при r = 0. (2.3.18) Условие на бесконечности имеет вид:

r2 + z 2.

A 0 при (2.3.19) Функции A, E, связаны соотношением i0 E + A + = 0. (2.3.20) c r Функции Hr, Hz выражаются через A следующим образом:

r = A, Hz = 1 (rA ).

H (2.3.21) z r r Из закона Ома следует, что = E.

j (2.3.22) Обозначим D0 образец с границей, проводимость которого будем считать конечной и равной его проводимости в жидкой фазе m, Di бесконечно проводящий индуктор с границей i, Db воздушное про странство между образцом и индуктором с проводимостью, равной нулю.

Введем функции F0 (r, z), Fb (r, z), Fi (r, z) следующим образом:

F0 (r, z), при (r, z) D0, A = (2.3.23) F (r, z), при (r, z) Db, b F (r, z), при (r, z) D.

i i Согласно [53] E = 0, / = i0 c1 при (r, z) Di, (2.3.24) где c1 = const. Тогда из (2.3.20), (2.3.23), (2.3.24) получим cc1 c Fi = =, (r, z) Di. (2.3.25) r r Отсюда и из 1-го условия (2.3.17) получаем условие c Fb |i =, (2.3.26) r которое можно рассматривать как граничное для функции Fb. Константа c должна быть выбрана из условия, что мощность выделяющегося в образце в единицу времени джоулева тепла равна мощности индуктора.

Вне индуктора отсутствует внешнее напряжение, поэтому [53] 0, (r, z) Di.

/ (2.3.27) Кроме того, в Db = 0, поэтому 2 Fb 1 (rFb ) + = 0, (r, z) Db. (2.3.28) z r r r Считая толщину скин-слоя в образце m = c/ 2m 0 = 0.0285 см (2.3.29) малым параметром, получим, что асимптотически при m 0 функция Fb удовлетворяет условию:

Fb | = 0. (2.3.30) (Это следствие первой формулы (2.3.17) на и формул (2.3.32), (2.3.33).) Кроме того, r2 + z 2.

Fb 0 при (2.3.31) Задача (2.3.26), (2.3.28), (2.3.30), (2.3.31) служит для определения функции Fb.

В образце электромагнитное поле проникает на глубину порядка m (см.

например [37, с. 408]). Главный член асимптотики ограниченного при n решения уравнения (2.3.16) при m 0 есть F0 (s, n) = a(s)e(1+i)n/m, (2.3.32) где n внутренняя по отношению к образцу нормаль к, s вектор касательной к, получающийся поворотом n на 90 градусов по часовой стрелке. Функция a(s) находится из 2-го условия (2.3.17):

m Fb m a(s) = = Hs, (2.3.33) 1 + i n 1+i где Hs | касательная к компонента амплитуды колебаний во времени вектора H, взятая на кривой. Это действительная функция, так как действительной является функция Fb. При n e(1+i)n/m.

Hs (s, n) = Hs (2.3.34) Функция в образце находится из условий (2.3.22), (2.3.20) с учетом j (2.3.27), (2.3.32), (2.3.33):

(s, n) = (1 + i)m m 0 Hs e(1+i)n/m.

j (2.3.35) 2c Найдем Hn в D0 :

1 (rF0 ) 1 (ra(s)) (1+i)n/m Hn (s, n) = = e = r s r s m (1 i) (rHs | ) (1+i)n/m = e. (2.3.36) 2r s Рассчитаем Hs (s, n, t), Hn (s, n, t), j (s, n, t) :

Hs (s, n, t) = Re[Hs ei0 t ] = Hs en/m cos(0 t n/m ), (2.3.37) Hn (s, n, t) = Re[Hn ei0 t ] = m (rHs | ) n/m = e (cos(0 t n/m ) + sin(0 t n/m )), (2.3.38) 2r s j (s, n, t) = Re[ ei0 t ] = j m m Hs en/m (cos(0 t n/m ) sin(0 t n/m )).

= (2.3.39) 2c Объемная плотность средней по периоду колебаний тока во времени мощ ности выделяющегося в образце джоулева тепла имеет вид [19, с. 19] 2/ 0 (Hs | )2 2n/m j q (s, n) = dt = e. (2.3.40) 2 m Средняя по периоду колебаний тока во времени пондеромоторная сила имеет координаты fs, fn, где 2/ (Hs | )2 2n/m 0 j H s fn (s, n) = dt = e, (2.3.41) 2 c 8m 2/ 0 j H n fs (s, n) = dt 0. (2.3.42) 2 c Проинтегрировав (2.3.40) и (2.3.41) по n, получим выражения для интен сивностей источников джоулева тепла q и магнитного давления f на n единицу поверхности образца:

0 m (Hs | )2, q (s) = q (s, n)dn = (2.3.43) (Hs | )2.

fn (s) = fn (s, n)dn = (2.3.44) Мощность индуктора N0 определяется из решения одномерной задачи для температуры [13] N0 = 4(rs fbs 1 s Tm /5)1/2 + 2fbm 1 rm Tm 2l = 9.442 · 1010 эрг/с, (2.3.45) 3 5 где rs = 5 см радиус образца, rm = 3 см средний радиус плавающей зоны. Эта мощность должна равняться мощности всего выделяющегося в образце джоулева тепла:

N0 = 2r(s)q (s)ds. (2.3.46) При решении задачи об определении параметров электромагнитного поля в полной постановке равенство (2.3.46) должно регламентировать значение константы c1 в (2.3.26). Однако мы не будем решать задачу для функции Fb, а ограничимся заданием функции f (x) на свободной границе плаваю щей зоны m, определенной формулой (Hs | )2 = H0 f (s/sm ), (2.3.47) где H0 максимум функции Hs | на кривой m, sm безразмерная длина m, s безразмерная длина дуги на m, отсчитываемая от нижней точки m (нижней точки трехфазного контакта). Из определения f (x) следует, что ее максимум равен 1, а областью определения служит отрезок [0, 1].

Предполагая, что мощность токов, выделяющаяся на m, равна 0.6N0, на основании (2.3.46), (2.3.43), (2.3.47) получим:


sm H0 = (8 · 0.6N0 /(l2 0 m r(s)f (s/sm )ds))1/2, (2.3.48) где r(s) безразмерная координата r соответствующей точки на m.

Перепишем формулы (2.3.40), (2.3.41) в терминах функции f (x) в без размерном виде:

0 H q (s, n) = fn (s, n), (2.3.49) (s, n) = H0 fn (s, n), fn (2.3.50) 8m где fn (s, n) = f (s/sm )e2n/Em. (2.3.51) В дальнейшем при расчетах будем рассматривать 2 варианта задания функ ции f (x). При вычислении формы m в разделе 1.2 эта функция задава лась в виде:

1/ 1 cos f (x) = ·54 · 4 3, ch(0 (0, 75x 0, 5)) cos где 0 = 0.1, 0 = 9 град. Она принимается в качестве 1 варианта задания f (x). Ее график кривая 1 на рис. 5.

Рис. На рис. 2 работы [54] для радиуса кристалла 3.8 см приведены получен ные расчетом форма свободной границы плавающей зоны и распределение по ней плотности электрического тока J. Отсекая верхнюю часть свобод ной границы (чтобы получить ее форму близкой к приведенной в [56]) и учитывая, что f пропорциональна J, получим функцию f (x), приве денную в табл. 3. Ей соответствует ломаная 2 на рис. 5.

Чтобы получить непрерывную функцию, разложим эту таблично задан Таблица k 0 1 2 3 4 5 6 7 f · 100 4.0 5.3 9.1 22.9 35.0 61.0 84.9 97.0 x · 100 0.0 1.7 6.3 13.7 18.0 23.0 29.4 36.0 44. k 9 10 11 12 13 14 15 16 f · 100 98.7 92.7 85.0 68.0 35.3 21.2 9.3 6.4 3. x · 100 52.5 60.9 67.2 76.1 85.3 87.8 92.4 96.2 ную функцию по полиномам Лежандра Xn (x) [12], удовлетворяющим усло вию нормировки Xn (x)dx = 1, n = 0, Nx.

Для этого определим коэффициенты bn по формулам:

(xk+1 xk ) bn = (fk Xn (xk ) + fk+1 Xn (xk+1 )), n = 0, Nx, k= где k номер колонки в табл. 3. Заметим, что нельзя брать число Nx боль шим из-за малого числа разбиений отрезка 0 x 1. Расчет показывает, что при Nx = 7 кривая Nx f (x) = bn Xn (x) n= хорошо аппроксимирует таблично заданную функцию внутри отрезка x 1, но вблизи его границ появляются участки немонотонности f (x).

Чтобы их ликвидировать, сделаем преобразование переменной x в правой части предыдущего равенства:

Nx f (x) = bn Xn (0.97(x + 0.05)/1.05).

n= Сделаем также растяжение вдоль оси ординат, обеспечивающее максимум f (x), равный единице. Так определенная функция f (x) принимается в качестве 2-го варианта задания пондеромоторной силы на m. Ее график кривая 3 на рис. 5.

2.4. Гидродинамическая часть задачи Пусть функции r(x, y), z(x, y) осуществляют конформное отображение прямоугольника 0 x 1, 0 y Y на область, занятую распла вом, так, что стороны x = 0, x = 1, y = 0, y = Y прямоугольника переходят соответственно в ось симметрии (r = 0) 0, свободную (пра вую) границу m, фронт кристаллизации c (нижнюю границу расплав – монокристалл) и фронт плавления f (верхнюю границу расплав – заготов ка). Пусть u(x, y, t), v(x, y, t), w(x, y, t) компоненты скорости расплава в направлениях x, y,. Введем функцию тока и вихрь по формулам 1 1 1 (uH) (vH) u=, v=, = 2, (2.4.1) rH y rH x H y x где (r/x)2 + (r/y) H= коэффициент Ламэ.

Форма уравнений осесимметричного движения жидкости с переменной вязкостью следует из тождеств 1 1 w w r3 r (divP) = 2 2 + = m rH x x r y y r 1 r r = 2 rw 2 rw + r2H 2 x x y y + r (rw) + r (rw), x x y y 1 (H(divP)x ) (H(divP)y ) = m H 2 y x 1 1 1 A B = (r) + (r) + + = H 2 x r x y r y x y 1 r r = 2 +r + 2 +r + H 2 x x x r y y y r A B + r + r + +, x x r y y r x y где 2 u 2 H 2 v 2 H A= +2 v +2 u, x H y H x y H y H x 2 u 2 H 2 v 2 H B= +2 v +2 u, x H x H y y H x H y P тензор напряжений. Эти тождества устанавливаются прямым вычисле нием. При этом используются формулы для компонент P согласно [18] и компонент divP согласно [22]. Формулы проверены в пакете "Математика" для конкретного случая, когда r = ex cos y, z = ex sin y, H = ex.

Введем модифицированные скорости U, V, W и модифицированный вихрь по формулам W = rw, =, (2.4.2) r U=, V =. (2.4.3) y x Поставим задачу для определения поля скоростей и температуры в без размерных переменных. Для должно выполняться уравнение импульса 1 1 r + 2 +r +U + rH 2 x t Re x x 1 r + 2 +r +V y Re y y 1 = G1, r + r (2.4.4) Re x x y y где r W 2 r W 1 G= + + rH 2 y r3 x r x y 2 1 1 1 H + rH 3 x x Re x x H y rH y 1 1 1 H + + rH 3 x y y H y rH x 1 1 1 H + rH 3 y x y x H x rH y 1 1 1 H + rH 3 y y y H x rH x Gr r r T + T Eum (Hfnxy ), (2.4.5) Re2 x x y y y y fnxy (x, y) = fn (s(y), n(x, y)), s(y) = H(1, y )d, n(x, y) = H(1, y)(1 x), y fn (s, n) определена в (2.3.51).

Выведем граничные условия для. На оси симметрии должно выпол няться условие симметрии = 0, x = 0. (2.4.6) x На свободной границе m условие Марангони [1, с. 195] · P · =, s n s где внешняя нормаль к границе, вектор, получающийся поворо n s том на 900 по часовой стрелке. Или n 1 T T pxy = =. (2.4.7) H y H y Согласно [18] 1 u 1 v u H v H pxy = m + 2 2. (2.4.8) H y H x H y H x Используя (2.4.1), (2.4.2) и условия u = u/y = 0 на m, представим pxy в виде 2v H pxy = m r +, x = 1.

H 2 x Подставляя это выражение в (2.4.7), обезразмеривая и пренебрегая членом, содержащим множитель v, ввиду большого отношения параметров M a и Re, окончательно получим:

Ma 1 T = ·, x = 1. (2.4.9) Re rH y Член, содержащий множителем v, мы отбросили для того, чтобы в даль нейшем точно расщепить задачи для вычисления и.

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

· P · = p = k · ( ), y = 0, y = Y, s n s v u xy где скорость стенки (твердого кремния), v u скорость расплава, k коэффициент проскальзывания. Используя (2.4.1), (2.4.2), (2.4.8), предста вим pxy в виде 2 v 2u H pxy = m r + 2.

H x H y Компонента v скорости на c,f равна проекции скорости стенки на u направление y, умноженной на S 1 z v = vc,f S ·.

H y Используя определение H и условия Коши - Римана, находим, что v z 1 H = vc,f S Kc,f, = Kc,f, H 2 y x x c,f где 2 z r 2 r z Kc,f =3 (2.4.10) x2 x x2 x H c,f кривизна границ фазового перехода. Кроме того, · = i vc,f z, · = i u, sv su c,f c,f H x где ic = 1, if = 1. Проведя обезразмеривание, получим условия про скальзывания в следующем виде:

1 Vc z = (2Kc Al) + (2Kc S Al), y = 0, r2 H y rH x 1 Vf z = (2Kf + Al) + (2Kf S + Al), y = Y, (2.4.11) r2 H y rH x где Kc,f определено в (2.4.10), Al безразмерный параметр проскальзы вания, kl Al =.

m m (В дальнейшем при тестовых расчетах полагается Al = 100, при расчете конвекции Al = 5000.) Здесь учтено, что на c,f безразмерная вязкость = 1.

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

Можно показать, что главные при h 0 члены в этих условиях определя ют условия проскальзывания с непостоянным Al = 1/(hH), где h шаг сетки в направлении, ортогональном к рассматриваемой границе.

Задача для функции тока 1 = rH 2, + (2.4.12) x r x y r y Rf = 0, x = 0, = Vf S, x = 1, r r = Vc S, y = 0, = Vf S, y = Y. (2.4.13) 2 Функция W должна удовлетворять уравнению импульса W 1 2 r 2 r + +U W + +V W rH 2 x t Re x y Re y 1 W W r + r =0 (2.4.14) Re x x y y и граничным условиям:

W = 0, x = 0 (2.4.15) это следствие из (2.4.2), условию отсутствия касательного напряжения в направлении на свободной границе 1 (W/r) W r px = 2 = H x r H x [18]. Или W r W 2 = 0, x = 1 (2.4.16) x x r и условиям прилипания на твердых границах W = c r2, y = 0, W = f r2, y = Y. (2.4.17) Для T должно выполняться уравнение энергии T 1 (U T ) (V T ) 1 T T + + r + r = Qefnxy rH t x y Pe x x y y (2.4.18) и граничные условия: условие симметрии T = 0, x = 0, (2.4.19) x закон Стефана - Больцмана на свободной границе 1 T T + Bi = 0, x = 1 (2.4.20) H x T и условия первого рода на фронтах плавления и кристаллизации T = T0, y = 0, y = Y. (2.4.21) Зададим также начальные условия = b, = b, W = Wb, T = Tb, t = 0. (2.4.22) Итак, задача (2.4.3), (2.4.4), (2.4.6), (2.4.9), (2.4.11) – (2.4.22) служит для определения неизвестных функций,, U, V, W, T. Компоненты скоро сти u, v, w восстанавливаются с помощью первых двух равенств (2.4.1) и первого равенства (2.4.2).

2.5. Численный алгоритм Введем обозначения 1 r 1 r F 1 =, p1 =, q1 = 2 +r 2 +r, Re x x Re y y r 2 r 2 r r µ1 =, F 2 = W, p2 =, q 2 =, µ2 =, Re Re x Re y Re r G2 = 0, F 3 = T, p3 = 0, q 3 = 0, µ3 =, G3 = Qefnxy. (2.5.1) Pe Тогда уравнения (2.4.4), (2.4.14), (2.4.18) можно записать в универсальной форме + Li + Li )F i = Gi, i = 1, 2, 3, ( (2.5.2) 1 t где 1 ((pi + U )F i ) i i F Li F i = µ, rH 2 x x x 1 ((q i + V )F i ) i i F ii L2 F = µ. (2.5.3) rH 2 y y y Пусть задана неравномерная сетка 0 = x0 x1... xN 1 xN = 1, 0 = y0 y1... yM 1 yM = Y.

Обозначим xn1/2 = (xn1 + xn )/2, hxn1/2 = (xn xn1 ), n = 1, N, ym1/2 = (ym1 + ym )/2, hym1/2 = (ym ym1 ), m = 1, M, hxn = (xn+1 xn1 )/2, n = 1, N 1, hym = (ym+1 ym1 )/2, m = 1, M 1.

Будем считать, что hxn+1/2 hxn1/2 = O(h2 xn+1/2 ), hym+1/2 hym1/2 = O(hym+1/2 ).

Введем также равномерную сетку на полуоси t 0 с шагом. Примем для любой функции f i (x, y, t) обозначение f i (xn, ym, k ) = fnm, ik причем индексы n, m могут быть дробными, а индексы n, m, i или k отсутствовать, если f не зависит от x, y, i или t или когда они несущественны.

В работе [7] для дифференциального оператора F µ F LF = U x x x получено разностное представление n = Un1/2 (Fn Fn1 ) + Un+1/2 (Fn+1 Fn ) In, F 2 hxn hxn где µn+1/2 (Fn+1 Fn ) In = n+1/2 cth n+1/2 hxn+1/ hxn µn1/2 (Fn Fn1 ) n1/2 cth n1/2, (2.5.4) hxn1/ hx U n1/2 =.

2µ n1/ Рассмотрим дифференциальное тождество (U F ) F + F U.

=U x x x Аналогично ему имеет место разностное тождество 1 (Fn+1 + Fn ) (Fn + Fn1 ) Un+1/2 Un1/2 = 2 hxn Un+1/2 (Fn+1 Fn ) Un1/2 (Fn Fn1 ) (Un+1/2 Un1/2 ) = + + Fn.


2 hxn hxn hxn Исходя из этих двух тождеств для дифференциального оператора (U F ) F LF = µ x x x получим следующее разностное представление:

1 (Fn+1 + Fn ) (Fn + Fn1 ) Fn = Un+1/2 Un1/2 In, 2 hxn где In определено в (2.5.4).

С учетом этого аппроксимируем дифференциальные операторы Li, Li 1 разностными операторами i, i следующим образом 1 ik Fnm = ik k ik ik xn [(Un+1/2m + pn+1/2m )Fn+1/2m 1 2h rnm Hnm µik ik ik n+1/2m n+1/2m cth n+1/2m k pik ik ik ik (Un1/2m + n1/2m )Fn1/2m (Fn+1m Fnm )+ hxn+1/ µik ik ik n1/2m n1/2m cth n1/2m ik ik + (Fnm Fn1m )], hxn1/ n = 1, N 1, m = 1, M 1, i = 1, 2, 3, ik Fnm = ik k ik ik ym [(Vnm+1/2 + qnm+1/2 )Fnm+1/ 2 2h rnm Hnm µik ik ik nm+1/2 nm+1/2 cth nm+1/ k ik ik ik ik (Vnm1/2 + qnm1/2 )Fnm1/2 (Fnm+1 Fnm )+ hym+1/ µik ik ik nm1/2 nm1/2 cth nm1/2 ik ik + (Fnm Fnm1 )], hym1/ n = 1, ni, m = 1, M 1, i = 1, 2, 3, где n+1/2m = hxn+1/2 (Un+1/2m + pik ik k ik n+1/2m )/(2µn+1/2m ), nm+1/2 = hym+1/2 (Vnm+1/2 + qnm+1/2 )/(2µik ik k ik nm+1/2 ), (2.5.5) n1 = N 1, n2 = n3 = N.

1 1 Распространим операторы Li, Li и 1, i при i = 1, 3, когда при i 1 2 x = 0 ставятся условия второго рода (2.4.6), (2.4.19), и функцию G1 на случай x = 0. Переходя к пределу (по правилу Лопиталя) при x 0, в i (2.5.3) с учетом µ = 0, U = 0, Fx = 0, H = rx при x = 0 получим (нижний индекс обозначает взятие производной по соответствующей переменной) Li F i = [((pi + Ux )F i )x ((2µi pi )Fx )x ], i 1 x x rx Li F i = [((qx + Vx )F i )y (µi Fy )y ], x = 0, i = 1, 3.

i i 2 x rx Продолжим сетку xn на один узел влево и учтем, что функции px, Ux нечетные, а функции µx, p, F, rx четные по x. Получим следующую разностную аппроксимацию оператора L1 :

(pik k ik (2µik ik ik ik x1/2m + Ux1/2m )F1/2m x1/2m p1/2m )(F1m F0m ) ik F0m ik =.

1 3 h rx0m hx1/2 x1/ (2.5.6) Аппроксимация оператора L2 при x = 0 имеет вид:

x ik F0m ik k ik ik = 3 [(Vx0m+1/2 + qx0m+1/2 )F0m+1/ r1m hym µik ik ik x0m+1/2 1m+1/2 cth 1m+1/ k ik ik ik ik (Vx0m1/2 + qx0m1/2 )F0m1/2 (F0m+1 F0m )+ hym+1/ µik ik ik x0m1/2 1m1/2 cth 1m1/2 ik ik + (F0m F0m1 )], m = 1, M 1, i = 1, 3.

hym1/ Функция G1 при x = 0 имеет следующий предел 1 2 Gr G1 = [I1 + I2 + (I3 I4 + I5 I6 + I7 I8 I9 I10 ) 2 (I11 + I12 )], rx Re Re где 2 Wxx ryx 1 Wxx xx yxx 2xx rxxx xx I1 =, I2 =, I3 =, I4 =, 3 2 2 2rx 4 rx rx rx rx y y yxx y rxxx xx y 1 xxxx 4 xx rxxx I5 = +, 2 2 2 rx rx rx rx 3 rx 3 rx y y y rxxx yxx 1 xx yxx xx rxy xx I6 =, I7 =, I8 =, 4 3 rx 2 rx rx y y y 1 xxxx 4 xx rxxx 1 y rxy yxx I9 =, I10 =, 2 3 rx 3 rx 3 rx 2 rx y y I11 = zy Txx + zyxx T ;

I12 = [zxx T ]y.

Причем с точностью O(x2 ) r1m 6r1m 6r2m rx =, rxxx = +, x1 (x2 x2 ) x2 (x2 x2 ) x1 2 1 2 21m 241m 242m 2W1m xx =, xxxx = 2 2 +22, Wxx =, x2 x1 (x2 x2 ) x2 (x2 x2 ) x 1 1 1 2(1m 0m ) 2(z1m z0m ) 2(T1m T0m ) xx =, zxx =, Txx =.

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

Граничные условия (2.4.16), (2.4.20) запишем в универсальной форме F i pi i (i 2)BiH i = iF F, x = 1, i = 2, 3.

x µ T С учетом этого распространим операторы i, i = 2, 3 на случай n = N [29, с. 45]:

ik FN m = ik k + pik 1/2m )FN 1/2m + ik [(UN 1/2m 1 N 2h rN m HN m xN 1/ ik ik ik (i 2)BiHN m µik m ik4 µN 1/2m N 1/2m cth N 1/2m ik ik N + FN m + (FN m FN 1m )], T0 hxN 1/ i = 2, 3.

Аппроксимируем (2.5.2) следующим образом ik+1 ik 2 2 ik ik (Fnm Fnm ) + (ik + ik )(Fnm Fnm )+ ik+1 ik (E + 1 2 ) +(ik + ik )Fnm = Gik, n = ni, ni, m = 1, M 1, i = 1, 2, 3, ik (2.5.7) 1 2 nm 0 тождественный оператор, ik оператор, отличающийся от ik где E 1 только при i = 3, n = N :

ik ik ik N m = k + pik 1/2m )N 1/2m + [(UN 1/2m 1 N 2h rN m HN m xN 1/ µik 1/2m N 1/2m cth N 1/2m ik ik 4(i 2)BiHN m µik m FN m ik ik N N + N m + T0 hxN 1/ ik ik (N m N 1m )], i = 3, весовой коэффициент, n1 = n3 = 0, n2 = 1.

[0, 1] 0 0 Уравнение (2.5.7) можно переписать в факторизованном виде:

ik+1 ik ik )(E + ik ) (Fnm Fnm ) = Gik (E + 1 2 nm (ik + ik )Fnm, n = ni, ni, m = 1, M 1, i = 1, 2, 3.

ik 1 2 0 Поэтому переход с k-го на k + 1-й временной слой при фиксированном i может быть осуществлен по алгоритму:

1) Определяем невязку nm = Gik (ik + ik )Fnm, n = ni, ni, m = 1, M 1.

ik ik nm 1 2 0 ik 2) Для всех m = 1, M 1 решаем системы уравнений относительно nm ik 0m = 0, при i = 2, 1 ik (E + ik )nm = nm, n = ni, ni, ik 0 ik+1 ik ik = (E + ik ) (g1m g1m ), при i = 1, N m 1k где g1m аппроксимация правой части (2.4.9).

3) При всех n = ni, ni решаем системы уравнений относительно поправ 0 ik ки nm ik+1 ik (g2n g2n ) ik n0=, (2.5.8) ik (E + ik )nm = nm, m = 1, M 1, ik ik+1 ik (g3n g3n ) ik nM =, (2.5.9) ik ik где g2n, g3n аппроксимация правых частей (2.4.11) при i = 1, (2.4.17) при i = 2 и (2.4.21) при i = 3. Далее полагаем ik+1 ik (g1m g1m ) ik ik 0m = 0, если i = 2, N m =, если i = 1, m = 0, M.

4) Для всех n = 0, N, m = 0, M вычисляем ik+1 ik ik Fnm = Fnm + nm.

Описанная разностная схема является монотонной при достаточно малом и консервативной (см. [7] и гл. 3). Если функции U, V и µi одного поряд ка, то она превращается в схему с центральными разностями для конвектив ных членов, которая имеет второй порядок аппроксимации по пространству.

µi, то она близка к схеме для уравнений Эйлера с разностями Если U, V против потока и имеет первый порядок аппроксимации по пространству. По времени описанная схема имеет первый порядок аппроксимации.

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

Представим в виде (несколько иначе, чем в [10, 26]) k+1 = k+1 + Pnm (g2n g2n ) + Qk+1 (g3n g3n ), k+1 1k+1 1k 1k+1 1k (2.5.10) nm nm nm где функции Pnm, Qk+1 удовлетворяют условиям k+ nm (E + 1k )Pnm = 0, Pn0 = 1, PnM = 0, k+1 k+1 k+ (E + 1k )Qk+1 = 0, Qk+1 = 0, Qk+1 = 1.

2 nm n0 nM 1k Тогда при реализации алгоритма для следует положить Fnm = k, Fnm = k+1, а правые части формул (2.5.8), (2.5.9) при i = 1k+ nm nm заменить нулями.

Разностная задача для выглядит следующим образом:

ax k+1 + bx k+1 + ay nm1 + by k+1 cnm k+1 + k+ nm n1m nm n+1m nm nm nm+1 nm +dyk+1 k+1 + eyk+1 k+11 = fnm, n = 1, N 1, m = 1, M 1, (2.5.11) k+ nm nm n1 nM Rc k+1 k+ = 0, = Vc S, m = 0, M, 0m Nm r rn k+1 = Vc S, k+1 = nM Vf S, n = 0, N, (2.5.12) n0 nM 2 где 1 ax =, bx =, nm rn1/2m hxn1/2 hxn nm rn+1/2m hxn+1/2 hxn 1 ay =, by =, nm nm rnm1/2 hym1/2 hym rnm+1/2 hym+1/2 hym cnm = ax + bx + ay + by, nm nm nm nm k+1 Qk+1 Hnm rnm (Al + 2Kf n ) Pnm Hnm rnm (Al 2Kcn ) yk+1 nm dyk+1 =, enm =, nm 2 Hn0 rn0 hy1/2 HnM rnM hyM 1/ nm fnm = rnm Hnm (k+1 Pnm g2n Qk+1 g3n ) dyk+1 k+1 eyk+1 k+ k+1 2 k+1 1k 1k nm nm nm n0 nM 2 k+ rnm Hnm Vc Pnm (2Kcn S Al)(zn+10 zn10 ) 2hxn rn0 Hn rnm Hnm Vf Qk+1 (2Kf n S + Al)(zn+1M zn1M ) nm.

2hxn rnM HnM Для решения задачи (2.5.11), (2.5.12) используем подход, предложенный В. Г. Зверевым в работе [15]. Сначала рассмотрим так называемый полиней ный метод решения данной задачи (для простоты записи опустим верхний индекс k + 1):

l+1/2 l+1/2 l+1/ ax n1m cnm nm + bx n+1m = ay nm1 by l l+1/ nm nm+ nm nm nm l+1/ dy n1 ey l 1 fnm, n = 1, N 1, m = 1, M 1, (2.5.13) nm nm nM ay l+1 cnm l+1 + by l+1 + dy l+1 + ey l+1 1 = nm nm1 nm nm nm+1 nm n1 nm nM l+1/ = ax l+1 bx n+1m fnm, m = 1, M 1, n = 1, N 1. (2.5.14) nm n1m nm Недостатком этого алгоритма является то, что члены nm+1, nM 1 в формулах (2.5.13) и член n+1m в формулах (2.5.14) берутся с нижнего итерационного слоя, что сильно ухудшает сходимость алгоритма. Поэтому, следуя [15], предположим, что nm+1 = nm nm + nm n1 + nm nM 1 + nm, (2.5.15) n+1m = nm nm + nm.

(2.5.16) Подставляя эти выражения в (2.5.13), (2.5.14) и учитывая неявно nm, n1, получим l+1/2 l+1/2 l+1/ ax n1m (cnm by nm )l+1/2 + bx n+1m = ay nm nm nm nm nm nm l l+1/ (dy + by nm )n1 (ey + by nm )l 1 (fnm + by nm ), (2.5.17) nm nm nm nm nM nm ay l+1 (cnm bx nm )l+1 + by l+1 + dy l+1 + nm nm1 nm nm nm nm+1 nm n +ey l+1 1 = ax l+1 (fnm + bx nm ).

l+1/ nm (2.5.18) nm nM nm n1m Следуя [15], аналогично методу неполной факторизации, запишем урав нение (2.5.11) в двух вариантах:

ay nm1 (cnm (ax + bx ))nm + by nm+1 + dy n1 + ey nM 1 = nm nm nm nm nm nm = ax n1m bx n+1m + (ax + bx )nm fnm, (2.5.19) nm nm nm nm ax n1m (cnm (ay + by ))nm + bx n+1m = ay nm nm nm nm nm nm by nm+1 + (ay + by )nm dy n1 ey nM 1 fnm, (2.5.20) nm nm nm nm nm где коэффициент нижней релаксации, 0 1.

Подставляя в (2.5.19), (2.5.20) выражения (2.5.15), (2.5.16), получим ре курентные формулы для определения коэффициентов связи:

l nM 1 = 0, nM 1 = 0, nM 1 = 0, nM 1 = nM, nM 2 = kn ay 1, nM 2 = kn dy 1, y y nM nM l l nM 2 = 1 + kn (cnM 1 + ey 1 ), nM 2 = kn (by 1 nM 1 + fnM 1 ), y y l nM nM y ay y y y nm1 = bnm nm y+ dnm, nm1 = bnm nm y+ enm, nm nm1 =, cnm by nm nm cnm bnm nm cnm bnm nm by l + f l l nm1 = nm nm y nm, m = M 2, 2, n = 1, N 1, (2.5.21) cnm bnm nm где 1 x y, fnm = ax l l x l x l kn = nm n1m + bnm n+1m (anm + bnm )nm + fnm, cnM cnm = cnm (ax + bnm );

x (2.5.22) nm x N 1m = 0, l+1/2 = N m, N 2m = aN 1m, N 1m cN 1m l+1/2 l+1/ bx 1m N 1m + fN 1m ax l+1/2 N nm N 2m =, n1m =, cN 1m cnm bx nm nm l+1/2 l+1/ bx nm nm + fnm l+1/ n1m =, n = N 2, 2, m = 1, M 1, (2.5.23) bx nm cnm nm где l+1/2 l+1/ l+1/ = ay nm1 + by nm+1 (ay + by )l+1/2 + fnm nm nm nm nm nm l+1/2 l+1/ +dy n1 + ey nM 1 + fnm, cnm = cnm (ay + by ).

nm nm nm nm y Заметим, что выбор коэффициента kn в формулах (2.5.21) находит ся в нашем распоряжении. В формулах (2.5.22) он определен так, чтобы nM 2, nM 2, nM 2 и nM 2 вычислялись по тем - же формулам, что и nm1, nm1, nm1 и nm1 при m M 1.

Запишем алгоритм перехода с l - го на l + 1 - й итерационный слой при решении системы (2.5.11), (2.5.12), l = 1, 2,..., L.

1) Если l = 1, то по формулам (2.5.21) определяются коэффициенты nm, nm, nm, по формулам (2.5.23) коэффициенты nm.

2) По формулам (2.5.21) определяются l. nm 3) Для всех m = 1, M 1 решаются системы уравнений (2.5.17) с учетом граничных условий (2.5.12).

l+1/ 4) По формулам (2.5.23) определяются nm.

5) Для всех n = 1, N 1 решаются системы уравнений (2.5.18) с учетом граничных условий (2.5.12) [34, с. 90].

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

1k+1 1k+ После решения задачи для k+1 вычисляются функции g2n и g3n nm и по формуле (2.5.10) восстанавливается функция k+1.

nm 2.6. Тестовые расчеты Если некоторая дифференциальная задача имеет стационарное решение и аналогичным свойством обладает ее дискретный аналог, то можно иссле довать дискретное решение на сходимость к решению дифференциальной задачи. Если 1 и 2 нормы погрешности какой-либо из искомых функ ций дискретного решения на сетках N M и 2N 2M соответственно, то порядок сходимости для этой функции определяется формулой = log2 (1 /2 ).

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

Рассмотрим следующую задачу 2 +2 = G, x y x x y Re y 2 + =, при 0 x 1, 0 y 1, x2 y = = 0, y = 0 и y = 1, = 0, x = 0 и x = 1, y = 0, x = 0, = g(y), x = 1.

x Функции G(x, y) и g(y) подберем так, чтобы функция = T = sin2 x sin2 y являлась точным решением этой задачи. В дискретном ана логе вместо условий /y = 0 поставим условия Тома [10]:

2 (x, 0) = (x, y1 ), (x, 1) = (x, yM 1 ).

h2 h2 1/ y1/2 yM Эта задача решалась на последовательности сеток методом установления по времени с помощью разностной схемы, описанной в предыдущем разделе.

Для каждой из функций,, U = /y и V = /x вычислялись 2 нормы погрешности Ri (норма в C и в L2 ):

1 Ri dy)1/2, i = 1, 2.

i = max |Ri | и i = ( dx 0 При Re = 0.1 в обеих нормах имел место второй порядок сходимости, при Re = 100 первый. Причем в последнем случае он начинал проявляться при бльших числах разбиений, чем в первом. При Re = 10000 численное о решение не сходилось к точному, так как последнее становилось неустойчи вым, но наблюдалась сходимость "в себе" с порядками, лежащими в преде лах 0.6 1.0 для всех функций кроме. Для порядок, вычисленный по нормам i, был отрицательным, а порядок, вычисленный по нормам i, равен 0.4. Эти результаты получены при 4N = 4M = 160. Отличие порядков сходимости для некоторых функций от 1 при Re = 10000 может быть объяснено недостаточным числом разбиений и разрывностью функции в угловых точках.

Был также произведен расчет порядков сходимости "в себе" и, соответствующих нормам i и i погрешности в определении функ ций,, W, T, U, V при решении задачи, описанной в предыдущих разделах, но с некоторыми упрощениями. Так был произведен расчет при Re = 1, M a = 0, Gr = 0, Eum = 0, Qe = 0, c = 0, f = 0, вязкость переменная. Конформное отображение квадрата 0 x 1, 0 y 1 на Таблица 4N 4M Параметр Функция W T U V 80 80 1.76 1.94 1.95 1.82 1.65 2. 1.93 1.91 2.03 2.02 1.76 1. 80 160 1.90 2.32 1.87 1.99 2.01 1. 1.95 1.94 2.01 2.01 1.96 1. 160 160 1.62 1.90 2.01 1.91 1.57 1. 1.57 1.69 2.03 2.03 1.85 1. 160 320 1.84 1.78 1.94 2.01 1.99 1. 1.75 1.27 2.02 2.01 1.86 1. область, занимаемую расплавом, задавалось по формулам c0 sinh x c0 sin( 0.5 + y) r(x, y) =, z(x, y) =, cosh x cos( 0.5 + y) cosh x cos( 0.5 + y) c0 выбиралось из условия r(1, 0) = Rc. Тогда c H(x, y) =.

cosh x cos( 0.5 + y) Параметр проскальзывания Al был равен 100. Результаты приведены в табл. 4.

Из табл. 4 видно, что почти все порядки сходимости лежат в пределах 1.5... 2.3 и большинство близки к 2 (теоретическое значение). Исключение составляет порядок сходимости функции при 4N = 160, 4M = 320.

1k Возможной причиной является то, что коэффициент n+1/2m при n = близок к 2 (см. формулы (2.5.5), (2.5.1)), а при n = 0 равен нулю (см.

формулу (2.5.6)). Другой возможной причиной является влияние ошибок округления.

Был также произведен расчет порядков сходимости "в себе" при 4N = 4M = 120, Gr = 0, Eum = 0, c = 0, f = 0, вязкость постоянная, осталь ные параметры как в разделе 2.2, геометрия области и параметр проскаль зывания как в предыдущем примере, то есть расчет конвекции Марангони Таблица Параметр Функция T U V 0.21 0.15 1.02 0.28 0. 0.60 0.51 0.67 0.47 0. при реальном (то есть очень большом) числе Рейнольдса. Результаты при ведены в табл. 5.

Отличие большинства порядков от 1 (теоретическое значение) может быть объяснено такими причинами, как недостаточное число разбиений, 1k влияние ошибок округления и разрывность коэффициента n+1/2m при переходе n от 0 к 1.

2.7. Результаты расчета конвекции Конвекция расплава в плавающей зоне обусловлена такими факторами, как:

1) пондеромоторная сила (результат действия магнитного поля) 2) термокапиллярный эффект (зависимость коэффициента поверхност ного натяжения от температуры) 3) сила плавучести, обусловленная зависимостью плотности расплава от температуры 4) вращение образца 5) протягивание образца.

Расчет производился методом, описанным в разделе 2.5, на последова тельности ортогональных разностных сеток, построенных в главе 1. При расчете задавалось 6 различных вариантов исходных данных. В 1, 2, 3, 5 и вариантах функции пондеромоторной силы и объемной плотности источни ков джоулева тепла определялись в соответствии с 2-м вариантом задания функции f (x) (см. раздел 2.3), а в 4 варианте в соответствии с 1-м вариан том задания функции f (x). Протягивание образца и переменность вязкости учитывались во всех 6 вариантах расчета. 1-й, 2-й и 3-й варианты расчета соответствуют действию 1-го, 2-го и 3-го факторов, обуславливающих кон векцию, отдельно от других факторов. В 4-м и 5-м вариантах учитывалось действие всех факторов, а в 6-м варианте действие 3, 4 и 5-го факторов.

Опишем порядок расчетов. Пусть xO y правая система координат в плоскости независимых пространственных переменных задачи, такая, что функции r(x, y), z(x, y) конформно отображают прямоугольник : x 1, 0 y Y = 1.08732 на двумерную область, занятую расплавом.

Так что отрезки y = 0, y = Y, x = 0, x = 1 переходят соответственно в фронт кристаллизации c, фронт плавления f, ось симметрии 0 и свободную границу m.

Начальные приближения для функции тока и модифицированного вихря брались при расчете удовлетворяющими граничным условиям при y = 0, y = Y и линейными по y. (Для это обеспечивало и выпол нение граничных условий при x = 0 и x = 1, а для граничного условия при x = 0.) Начальное приближение для температуры T бралось удовлетворяющим граничным условиям при y = 0 и y = Y и параболи ческим по y с максимумом, равным T0 + 1 (это обеспечивает перегрев расплава 50K, граничное условие при x = 0 для T при этом также выполняется). Если решалась задача с вращением, то начальное приближе ние для модифицированной азимутальной компоненты скорости W бралось равным нулю, а угловые скорости вращения нижнего и верхнего цилиндров задавались зависящими от времени t по закону c,f t (t) = (1 et )c,f, t [0, 10).

При t 10 полагалось c,f t (t) c,f. Шаг по времени во всех расчетах был равен 103, что соответствует 1.5·104 с размерного времени, малость обусловлена малостью минимального шага сетки в переменных x, y : xN xN 1 = 4.275 · 103 /N, где N число разбиений по оси O x (см.

раздел 1.3), безразмерный параметр проскальзывания Al = 5000, весовой параметр = 0.99 (см. раздел 3.4), а коэффициент нижней релаксации при решении задачи для = 0.8.

Число итераций L при решении задачи для на каждом шаге по вре мени определялось следующим образом. Задавалось некоторое малое число 0. Производилось Lmin итераций. Если 2 H 2 dxdy/ f H 2 dxdy 0, 2 Hds + = n (r/x)2 + (r/y) где n функция невязки, H = коэффициент Ламэ, f функция правой части при решении задачи для, граница, s длина дуги, то итерации прекращались. В противном случае они продолжались до тех пор, пока не станет меньше 0 или число итераций не станет равным Lmax. Число итераций K по времени при решении всей задачи определялось аналогично числу итераций при решении задачи для с заменой Lmin, Lmax на Kmin, Kmax ;

0 на 0 и на (/t)2 H 2 dxdy/ 2 H 2 dxdy.

= Для всех вариантов было Lmin = 2, Lmax = 20, 0 = 0 = 108. Сначала задача для вариантов 1–3 решалась на сетке 30 30 с Kmin = 104, Kmax = 5 · 105, затем на сетках 30 60, 60 60 и 60 120, каждый раз с начальным приближением, построенным на предыдущей сетке, при Kmin = Kmax = 1.5 · 104. И окончательный расчет делался на сетке 120 120 при Kmin = 104, Kmax = 5 · 105.

Для вариантов 1, 2, 3 число итераций по времени составило на сетке 30 30 соответственно 76499, 199966 и 344737, а на сетке 120 120 67185, 358165 и 500000. Значение параметра в 3 варианте на сетке 120 составило 1.2 · 108. На последнем шаге по времени число итераций при решении задачи для у 1–3 вариантов на сетках 30 30 и 120 составляло L = 2. Значение начального приближения на последнем шаге по времени на сетке 30 30 для вариантов 1, 2, 3 составило 1.3 · 1011, 4.4 · 1012 и 3.2 · 1012, а на сетке 120 120 1.2 · 1011, 1.5 · и 2.9 · 1013 соответственно. После проведения 2 итераций значение составило на сетке 30 30 для вариантов 1, 2, 3 соответственно 9.9 · 1016, 5.6 · 1015 и 7.7 · 1015, а на сетке 120 120 3.0 · 1013, 1.0 · и 1.5 · 1013.



Pages:   || 2 |
 





 
© 2013 www.libed.ru - «Бесплатная библиотека научно-практических конференций»

Материалы этого сайта размещены для ознакомления, все права принадлежат их авторам.
Если Вы не согласны с тем, что Ваш материал размещён на этом сайте, пожалуйста, напишите нам, мы в течении 1-2 рабочих дней удалим его.