Содержание к диссертации
Введение
ГЛАВА 1. Построение ортогональной разностной сетки 11
1.1. Пос і роение одномерной сетки, сгущающееся на краях 01 резка М
1.2. Вычисление формы свободной границы 1С
1.3. Построение пеорюгоналыюй двумерной
ГЛАВА 2. Моделирование конвекции при бестигелыюй зонной плавке в магнитном поле 28
2.1. Размерные парачеірьі задачи 29
2 2. Безразмерные кршерии подобия 31
2.3. Элекіромагшпнам часіь задачи 32
2.-1. Гидродинамическая часів задачи 11
2 5. Численный алюритм 17
2.0. Тес і оные расчет 58
2.7. Результат рас чопі конвекции G1
ГЛАВА 3. Условия монотонности факторизованиой разностной схемы для параболического уравнения 78
3.1. Описание разностной схемы 79
3.2. Доспи очные условия моноюниосіи в общем случае 81
3.3. Необходимые и достаточные условия моноюниости для задачи с малым числом разбиений 85
3.4 Пример расчет 90
ГЛАВА 4. Задача о вращающемся контейнере 91
4 1. Постановка задачи 91
4.2. Аеимпюшка вихря и функции юка в окрсччноешючек контакта 102
4.3. Аналитическио решения 109
4.1 Меюд и результаты расчпа течения 115
Заключение 125
Список литературы
- Вычисление формы свободной границы
- Безразмерные кршерии подобия
- Необходимые и достаточные условия моноюниости для задачи с малым числом разбиений
- Аналитическио решения
Введение к работе
Актуальность темы. Зонная плавка — кристаллографический метод рафинирования материалов, состоящий в перемещении узкой расплавленной зоны вдоль длинного твердого стержня из рафинируемого материала. В диссертации рассматриваются два варианта зонной плавки.
1 вариант: бестигельная зонная плавка (БЗП) в магнитном поле (МП),
применяемая, в частности, для выращивания монокристаллов кремния
большого радиуса (более 5 см). Верхняя (заготовка) и нижняя (выращи
ваемый монокристалл) части вертикального цилиндрического образца
медленно движутся вниз и вращаются в противоположных направлени
ях с разными угловыми скоростями. Часть нижней границы заготовки
покрыта жидкой пленкой, остальная часть граничит с плавающей зо
ной, находящейся между заготовкой и монокристаллом, поддерживаемой
в жидком состоянии неподвижным источником высокочастотного элек
тромагнитного поля — индуктором, и удерживаемой между твердыми
частями образца силами поверхностного натяжения и магнитного давле
ния. Токи, наводимые индуктором, сосредоточены в тонком скин-слое,
примыкающем к свободной границе расплава. Они приводят к выделе
нию джоулева тепла и создают пондеромоторную силу, направленную
ортогонально свободной границе, экспоненциально убывающую при уда
лении от нее, и являющуюся одним из источников конвекции в расплаве.
Полученные методом БЗП в МП монокристаллы кремния большого радиуса используются в основном в двух направлениях:
в силовой электронике — создание тиристоров, силовых транзисторов, используемых в мощных силовых преобразователях;
при изготовлении высокоэффективных солнечных батарей.
2 вариант: зонная очистка поликристаллического полупроводниково
го материала (в частности, германия). Цилиндрический полый контейнер
наполнен полупроводниковым материалом и расположен под небольшим
углом к горизонту. Он вращается и одновременно совершает медленное
поступательное движение вдоль своей оси, при этом часть его нагревает
ся до высокой температуры. Материал полупроводника плавится. Обра
зуются фронты плавления и кристаллизации. Расплав не целиком запол
няет контейнер: имеется свободная поверхность. Предполагается также,
что расплав отделяет от стенки контейнера тонкий слой мелкодисперс
ной смазки.
Очищенный полупроводниковый германий применяется для создания транзисторов и солнечных батарей.
Осесимметричная нестационарная задача о БЗП в МП была численно решена в полной постановке немецким ученым А. А. Мюльбауэром с
четырьмя соавторами в 1995 г. Были рассчитаны характеристики электромагнитного поля, найдены форма границы плавающей зоны, поля скоростей и температуры. При этом, чтобы не слишком мельчить сетку вблизи свободной границы, пондеромоторная сила была снесена из уравнения импульса в граничное условие для вихря. Это удалось сделать в предположении о малости конвективных членов по сравнению с вязкими в уравнении импульса в скин-слое. В 1997 г. А. А. Мюльбауэр с двумя латвийскими учеными рассчитывает распределение примеси в растущем кристалле при БЗП в МП, а в 1999 г. группой немецких и латвийских ученых решается аналогичная задача уже в трехмерной постановке. В 2001 и 2005 гг. решаются осесимметричная и трехмерная задачи о БЗП в МП при наличии дополнительного низкочастотного индуктора, позволяющего получить дополнительное управление процессом.
Однако анализ показывает, что указанное предположение, сделанное во всех перечисленных работах, на практике не выполняется. Поэтому актуальной является разработка модели процесса БЗП в МП с учетом пондеромоторной силы не в граничном условии, а в уравнении для вихря (уравнении импульса).
Во втором варианте зонной плавки течение является существенно трехмерным и весьма сложным. Прямой его расчет — очень трудная задача. Поэтому для эффективного расчета гидродинамики расплава область течения предлагается разделить на ядро, где продольная компонента скорости мала и движение в первом приближении можно считать плоским, и пограничные слои возле фронтов плавления и кристаллизации. В диссертации рассматривается только течение в ядре, то есть решается задача о плоскопараллельном стационарном движении вязкой несжимаемой жидкости в горизонтальном цилиндрическом вращающемся контейнере. В имеющейся литературе рассматривался лишь случай, когда жидкость целиком покрывает стенки контейнера, то есть область течения является двусвязной, а на твердой границе жидкость прилипает к стенкам. Поэтому актуальной является задача о движении жидкости во вращающемся контейнере для случая, когда область течения является односвязной, что соответствует экспериментальным условиям во втором варианте зонной плавки, а на границе с твердой стенкой ставятся условия проскальзывания, моделирующие действие слоя мелкодисперсной смазки.
Цель работы: осуществить математическую формулировку задач, разработать численные алгоритмы и произвести численные расчеты движения расплава при БЗП в МП и во вращающемся контейнере.
Научная новизна. В работе впервые
решена гидродинамическая часть задачи о БЗП в МП с учетом пондеромоторной силы в уравнении для вихря;
выведены уравнения осесимметричного движения жидкости с переменной вязкостью для вихря и азимутальной компоненты скорости в дивергентной форме в ортогональных переменных;
реализована консервативная монотонная разностная схема при решении задачи о БЗП в МП на ортогональной сетке, полученной с помощью конформного отображения прямоугольника на область, занятую расплавом;
численно и аналитически решена задача о движении жидкости во вращающемся контейнере для случая, когда область течения является в первом приближении сегментом круга.
Достоверность численных результатов обосновывается тестовыми расчетами с проверкой на сходимость к точному решению задачи, если таковое имеется, и на сходимость "в себе" в противном случае.
Научная и практическая ценность. Работа представляет собой существенное продвижение на пути к построению полных двумерных моделей процессов зонной плавки, которые могут использоваться для оптимизации этих процессов.
Методы исследования. При численном решении задач расчета конвекции при обоих вариантах зонной плавки использовался конечно-разностный метод переменных направлений. При построении асимптотики вихря и функции тока в окрестности угловых точек области и построении аналитических решений при втором варианте зонной плавки использовался метод интегральных уравнений. При аналитическом решении задачи для определения формы свободной границы — метод функции Грина.
Апробация работы. Результаты работы докладывались на Всесоюзном семинаре по математическому моделированию процессов кристаллизации (Рига, 1989 г.), на десятой Зимней Школе по механике сплошных сред (Пермь, 1995 г.), на третьей Сибирской школе-семинаре "Математические проблемы механики сплошных сред" (Новосибирск, 1999 г.), на Всероссийской конференции "Теория и приложения задач со свободными границами" (Бийск, 2002 г.), на семинаре отдела прикладной гидродинамики Института гидродинамики им. М. А. Лаврентьева СО РАН и на объединенном семинаре "Информационно-вычислительные технологии" под руководством академика Ю. И. Шокина и профессора В. М. Ковени в Институте вычислительных технологий СО РАН. .
Публикации. По теме диссертации опубликовано 7 печатных работ. Из них 4 — в изданиях, рекомендуемых ВАК для представления
Рис. 1
основных результатов диссертации (в журнале "Вычислительные технологии"), 1 — монография, в которой 5 глава написана лично автором, и 2 статьи в сборниках. Все публикации в журналах — без соавторов.
Структура и объем диссертации. Диссертация состоит из введения, 4 глав, заключения и списка литературы. Работа изложена на 135 страницах, содержит 19 рисунков. Список литературы содержит 67 наименований.
Вычисление формы свободной границы
Пусть А — длина дуги, оіечишваечая or нижней точки свободной границы Г„, плавающей зоны, # ,р — амплшуда колебаний во времени ь-fi компо-ношы напряженное ги МП на Г = Гш U Г;, Г; — граница жидкой пленки, q — поверхносіная плоіносіь псючников джоулева тепла; /„ — функция магии і HOI о давления, Лго 9.442-10ш эрг/с — мощное і ь индуктора, вычисленная для радиуса образца о см, WQ = 1.76 107 рад/с — круговая часчоїа іока в индукторе, єт = 2.85 0 2 см — юлщіша скип-слоя в жидкой фазе. И мою г мес і о соотношения (см. раздел 2.3): q{s) = ( ш/(1бтг))(/Дг)2, (1.2.1) /п(Я) = (1/(1б7г))(/Дг)2, (1.2.2) ЛГ0Й J 2irr(s)q(s)ds. (1.2.3) г Последнее равенсіво имееі место при условии, чю джоулевым іепловьіде лением в івердои )tue кремния можно пренебречь.
Введем функцию /{х), заданную на единичном о і ре же и имеющую единичный максимум, описывающую форму функций магпиіпого давления и джоулева тепловыделения: m = (Hi/(m))f(s/sm), (1.2,1) где Н{) — .максимум Я(гш, s — безразмерная длина дуги, sm — безразмерная длина кривой Гт. В качестве масштаба длины выберем I — 1.5 см. Предположим, что мощноеп токов, выделяющаяся на Г7„, равна O.GiVo Тогда на основании (1 2.1)-(1.2.4) для //(, получим следующее выражение. Щ = (8 0,6AV(/Wm J ФШ АтМ ))1 . (1.2.5) о Сисіема уравнений для функции r(s), z(s), парамеїрически задающих іраницу Г,„, и безразмерных переменных имеег вид [56. 59]: r (.s) = sin 0(A), z (s) = cos B(s), e (s) = + Bo(z(s) + Fef(s/sm) - со), (1.2.6) fl. .ft.i (1.2.7) іЬтгршдІ m где 4 2 pm = 2.53r/eu — ііло і нос гь жидкого кремния, g = 980см/е — ускоре-ішс свободного падения, ат = 737 дин/ем коэффициент поверхноепюю наїяження расплава.
К (1.2 С) нужно присоедини іь условия Копій: г(0) = Лг, z{Q) = Zet в{О) = 0п (1.2.8) где вс = 11 г])ад — уюл между касательной к свободной границе при 6 = 0 и направлением z, при коюром процесс БЗП проіекаеі усіойчиво (выращиваемый крисіалл имеет постоянный радиус [60]).
Для определения консіанг со, sm необходимо іюсіавить 2 дополни і ель-иых условия, например, фш) = Я/; z(sm) = Zf. (1.2.9) Функция f(x) должна, по видимому, ])езко возрастать от значения, близкою к нулю, в окрееіносіи точки х = 0, и заіем плавно меняіься до ючки х" = 1. Зададим эту функцию в виде /W= lch(/3b(0,75і-0,5))- 3 J 5"4J "4"3, где / = 0.1; іро — 9 ірад. Ее график показан на рис. 1. Расчеі формы свободной границы производился следующим образом Сначала задача (1.2,6) - (1.2.9) решалась при IIQ = 310.9 Гс (оценочное — 0.8 — 0.6 — 0.4 — 0.2 — 0 —j—! j j—у j f— 0 0.2 0.4 0.6 0.8 1 Рис. 1 значение Но). Зачем по формуле (1 2 5) было определено Щ = 287.7 Ге. После эюго был проишеден окончательный расчет формы свободной границы. Задача (1.2.G) - (1.2,9) решалась меюдом Эйлера с помощью визуальной пристрелки но парамеїру со- (Для каждого значения со с і роился график свободной границы, после чего значение CQ подправлялось по меюду деления огрсжа пополам.) Исходные данные были следующие. Rc = 3.4, Zc = 0, Єс = И град, Щ = 0.8884, Zf = 0.9020. Число раз-биений по s было принято равным 256. Резулыаш расчеюв следующие: Fe — 0.4428; Во = 7.569; sm = 3; с() = 0.801898. График свободной границы приведен в следующем разделе (рис. 3)
Безразмерные кршерии подобия
Задача ставится в переменных вихрь - функции іока в сооївен ішш с приближением Буссиниска. Дифференциальные уравнения для функции тока, вихря, азимуіальноіі компонент скоросіи и іемпера-іурьі записываю ієн в диверіеїшюй форме. Дифференциальные онера юры конвективною и диффузионного переноса аппроксимпруюіся моноюнны-\ш консервативными разностными операторами, имеющими вюрой порядок аппроксимации по пропранегвенным переменным п])и малых числах Рей-нольдса и Пекле и первый — при больших числах Ройнольдса и Пекле {здесь псиользуюіся результаты работы [7, с 280-288]). При решении разностных уравнений исіюльзуеіся меюд, позволяющий ючным образом разделить задачи вычисления вихря и функции кжа, описанный в рабо і ах 10, 26. При решении на каждом временном шаіе уравнения для функции і ока исполь-зуеіся эффективный итерационный алгоритм, не требующий информации о евойеівах и границах спектра разностного онера юра, обобщающий ал-горшм, описанный в рабо і е 15. Расчет производпіся для двух варианшв задания пондеромоторной силы на свободной границе и для них іюлучаюі-ея принципиально разные режимы конвекции. В первом вариаше максимум и минимум функции тока колеблю і ся с периодом 0.2 с и движение суще-сівенно нееіацнонарно Во вюром вариаше максимум и минимум функции тока изменяются со временем моноюнно и движение постепенно выходи і на режим, близкий к сіацпонарному. Основные резулыаш главы изложены в работах авюра [G2, Gl, Go, G7J.
Плопюеп» р и поверхноспюе натяжение и расплава будем считать линейными функциями температуры Т (перемешюеіь плоіносіи учишваеіся в приближении Буссинеска). Кроме того, р терпит скачек при фазовом переходе. Кинематическую вязкоеіь иредсіавим функцией вида v[T) = au + Kl{T-cu).
В іабл. 2 приведены названия, обозначения и значения материальных констант [32. 50]. Символ "т" означает жидкую, "я" —твердую фазу кремния. Коэффициент av вычисляемся по формуле (h = vm Kl rm -cv).
Пусіь Оrzip — цилиндрическая система координат, где г — полярный радиус. р — полярный угол. В сечении р = const ось Or проведем го-ризошалыю вправо через нижнюю точку трехфазного кон і акт а, а ось Oz — вер шкально вверх вдоль оси образца. Кроме перечисленных нам понадо-бяіся следующие размерные парамеїрьі: толщина скин-слоя єт = c/y/2njmWQ = 2.85 10"2см, радиус монокриеіалла ге = 5.1 см, координата г точки контакта расплава с жидкой пленкой (верхней правой угловой точки плавающей зоны) гj = 1,333 ем, скорость проіягивания монокристалла vc = —5.52-10"3 см/с, скороеіь ироіягивания заютовки vj — vc(rc/rf)2 = —8.085- Ю-2 см/с (здесь не учитываемся, что большая часть массы поступает в расплав из жидкой пленки, а величина Vf принята такой, какой она была-бы, если-бы заюювка имела радиус гу), угловая скорое і ь вращения верхней а у = —2.16 рад/с и нижней шс — 5.24- Ю-1 ])ад/с границ с[)азовою перехода [56], характерный размер и плавающей зоне (равный половине ее высоты) / = 1.5 см, харакіер-ная скорость движения расплава v\ — 10 см/с (см. раздел 2.7), характерный перегрев расплава AT = 50 К [56], мощность индукюра NQ = 9.442 10і" эрг/о (см. раздел 2.3).
Максимальная на Г,н амплиіуда колебаний «о времени касательной к Г,„ компонешы вектора напряженности магнитного поля при сделанных в разделе 2.3 предположениях определяемся по формуло Яо - {8 0.6AV(/W„. f r(s)f(s/sm)ds))V\ І) где s — безразмерная длина дуги на Г,„, оіечиїьшаемая от нижней іочки Гни Г(А) безразмерная координата г соогіюіегвующен точки Гш, sm — значение s в верхней точке Гт, f(x) — функция, описывающая форм) іюверхносшоіі пло і носі» исючнпков джоулеиа тепла на Гт, заданная на единичном о і ре же и имеющая единичный максимум.
Необходимые и достаточные условия моноюниости для задачи с малым числом разбиений
Здесь не выписаны члены с заведомо неоіршпиельньїчи коэффициент-ми п]ш F„±i,Ha-i. Пусіь 7 0- Тогда нужно ришиїь сискзму неравежчв оіносиїельно т : 1 - {1 - у)т(х + у) + ту 0, хе(0,Х}, у Є (О, У], г rain {(1-7)/(7 (1-т)/(72П}, где А , К определены в (3.2.1). Обозначим f(x,y) = 1-(1- J)T(X + 1/) + я/2т ху. Имеем: df/дх = —(1 — "/)т + 72т22/ 0 при у у0 = (1 - 7)/(72т) и 5// = -(1 - і)т + 72г2ж 0 при х т/о- Из 2-го нсравепсіва (3.2 4) имеем: Л" г/«; У г/ц. Следовательно, минимум функции f(x,y) досіигаеісн при х = X; y = Y. Теперь порвое пернвенсіво (3.2.-1) можно заменил, на 1-(1- 7)т(Х + Y) + 72г2ХУ 0. (3.2.5)
Ею решение есть т т12 (см. формулы (3.2.1)). Пусіь дли определенное! и X Y. Составим разность (г12 - T)2y2XY = (1 - 7)(А -Y) — \Z"f2{X — У)2 + {1 - 27)(А + У)2. Предположим, чю 7 чаково. чю подкоренное вы])ажение неоірицаїельно Тогда мы имеем разность двух неоірица іельньїх чисел, а последняя имеет чот - же знак, чю и ])а їіюсіь их квадрачов, коюрая раина (1 — 27)(2А 2 + 2У2). Таким образом, при 7 1/2 т12 г13, при 7 = 1/2 ти = т13, а при 7 1/2 г1 2 г13 на учасіке существования вещественною г12. Если же 7 чаково, чю г12 комплексное, чо неравенство (3.2 5) выполняемся при любом г и, следовательно, ограничение на г дает только функция т13. При этом "/ 1/2. Определим т1 при 7 — 0 как предельное значение т12 при 7 - 0- В резульчаїе получим (формулу для г11 в (3.2.1) (тіо несколько более грубая оценка чем ча, которая пепоередеівенно следует из (3.2.3)).
Расемоірим чіan 2 Мы имеем се]шіо одномерных задач: -Італії! + (1 + 7 JC1/2 - irh AZlin = І. n = W T, ,/2 = (Е + 7"\2)5i, 411/2 = № + 7гЛ2) /2и„ Ш = ТДТ Т.
И знеспю (см. [11, с 270]}, чю для монотонности схемы на г)іом мапе достаточно неоірицаїельносіи коэффициентов afim,b mi котрая имеем1 место is силу условий предложения, и свойства диагональної о преобладания-1 + утс т ут(а т + bfim), коюрое можно записать в виде г г21, где г21 вычисляется по формулам (3.2.1).
На этапе 3 опять имеем серию одномерных задач для коюрой аналогично зі any 2 получаем т т 22. Предложение докатно. Необходимые и достаточные условия монотонности для задами с малым числом разбиений
Доказанное в предыдущем разделе предложение даег ограничения на шаг г, выполнения которых достаточно для МОІЮЮННОСІИ схемы. Возникает вопрос: насколько зіи ограничения близки к необходимым и, в часгносш, нвляен я ли схема безусловно немонотонной при 7 — 1 {Ъплх — 0 при 7 — 1)-Для опима на ти вопросы расс.моірим разностную схему для уравнения іешюнроіюдносіи OF/Ot - X{02F/0x2 + d2F/dtf) = 0 с числами разбиений Лг = М = А и равномерными шагами по осям ж, у (здесь х — const — коэффициент температуронроводносіи). Положим на іранице обласчи F = 0. Исключив граничные тючки, как показано в [23, с. 3G3], получим разной ную задачу с 9-ю узлами на каждом временном слое
Аналитическио решения
Таким образом, решение задачи (4.1.4())-(4.1.42) даеіся формулами (4.4.15)-(4.4.21). При А1 оо функция Ф(х) имеет в точках контакіа ішіегрируемую особениоегь и, вследствие этого, /" тоже. Вмесіе с тем функции f(x) и f[x) в этих точках непрерывны. Для решения (4.3,33), (1.3 35), построенною при А1 = оо, особенноеіь функции Ф(х) в точках кошакіа не интегрируема (эю видно из (4.3.35)) Поэю.му шпеграл в правой часі и (4.4.15) расходится и решения задачи (4.1.4())-(4.1.42) не еущесінуег.
Приведем пример численного эксперимент при значениях пара-меіров г/о = 0.5 рад (еоотиеіеівенно h = 0.88, т. е. цилиндр заполнен почти полноеіью, /3 = 1.19 (см. (4.1.32), (4.2.2)), /1/ = 1.5. Парамеїрьі Re и Во будут варьироваться, а параметр Fr не рассматривается, так как функция f(x) зависит or него линейно. Обратим внимание на следующие особенное пі расчет.
1) Среди асимпюшческих условий наибольшую ошибку дает условие для вихря (4.2.38). Для достижения точноеіи 10 3 в ею выполнении следам браті» о Ю- 12.
2) После введения сеіки в облаеіи (, т?) Є [-)] х [іод] сооївеїсіву-ющая часть области D в плоское! и переменных (х,у) буде і разбит па ячейки, диаметры коюрых приближенно равны \/Д - Ч- Дг/2//( , г/), где Д, Дг/ — шаги сеюк по осям , г}. Функция Я имеет максимум при = 0, ц = ЇДЬ равный 3.9. Чтобы обеспечить малость максимума диамеїров ячеек, нужно либо делать очень большое число разбиении но оси (порядка 400, чюбы значение эшго максимума было меньше 0.2 при о — Ю), либо использовать неравномерную ееіку.
3) При вычислении шпегралов (4.4 15), (4 4,17) функции Ф, определенная формулой (1 1.13), должна бьиь доопределена вне о і режа fj о с помощью асимпюіических фо])мул (4.2.31). (4.2.36). При зі ом неизвесшые конпангы су1, с]1 приближенно определяюіся по решению внуі])и чю-ю о і резка. IГервый член в выражении (4.1.43) имеет порядок 0(е Щ, в юрой — 0(), треінй — 0(е(2 ІІ№). Ишегралы or первых двух членов по внешносш о і речка о пренебрежимо малы, а ишеграл or іреіьего члена — порядка единицы. При расчеіе функции f(x) он учніьшалея в виде авали іически вычисляемого слагаемого.
На риг. 17, а показаны линии тока, на рис. 17, б — изолинии вихря, носі роенные по решению задачи (4.1.33)-(4.1.37), рассчшапному при Re = 10. На риг. 18, а. б приведены аналогичные графики для Re = 150. Функции Ф и UJ неположительны. Максимумы их модулей, сооївеїпвемио, сосіав-ляюі 0.333 и 1.673 при Re = 10 и 0.335 и 1.671 при Re = 150. В ючках контакта скорое і ь жидкости раина нулю Максимум скоросіи в D сосіавляег 0.693 при Re = 10 и 0.764 при Re = 150. В первом случае он досіигаеіся на твердой стенке, во в юром — на свободной границе. Результаты получены при о = 10 и числах разбиений по осям , //, равных Л = 64, NTI = 32. Сегка при расчеіе была неравномерной. Так, шаг по оси увеличивался оі значения 0.03 при = 0 до значении с? 1.5 при = ц. Для уменьшении влияния схемной вязкое і и разносіи проіив поюка первого порядка при аппроксимации первых производных и были заменены разнос і ими прошв поюка вюрого порядка ючносш Все остальные члены в уравнениях и граничных условиях также аппроксимировались разпостыми формулами вюрою порядка, Сисіема решалась методом пяш-диаюнальной прогонки [31].
