CC BY
61ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ПРОЦЕССА ВЫРАЩИВАНИЯ КРИСТАЛЛА ОА:ОЕ ПО ВЕРТИКАЛЬНОМУ МЕТОДУ БРИДЖМЕНА В ПРИСУТСТВИИ ВРАЩАЮЩЕГОСЯ МАГНИТНОГО ПОЛЯ
Т.П. Любимова, О. А. Хлыбов
Институт механики сплошных сред УрО РАН,
614013, Пермь, ак. Королева, 1
Проведено численное моделирование процесса выращивания кристалла Ga:Ge по методу вертикального Бриджмена в двумерной нестационарной постановке с учетом искривления фронта кристаллизации. Исследовано влияние аксиально-симметричного вращающегося магнитного поля на форму фронта кристаллизации, структуру и интенсивность течения и сегрегацию примеси в расплаве. Получены нестационарные картины меридионального и азимутального течений при превышении критического значения величины магнитного поля, когда течение Куэтта теряет устойчивость. Показано формирование вихрей Тейлора.
ВВЕДЕНИЕ
В настоящее время основным промышленным методом выращивания полупроводниковых кристаллов является метод Чохральско-го (Cz) и его модификации (VCz и др.), которые позволяют получать кристаллы большого диаметра (до 300 мм). Однако в некоторых случаях более перспективными оказываются другие методы, такие как жидкая зона (FZ), движущийся нагреватель (Traveling heater technique), вертикальный метод Бриджмена (VB) и др. Последний в основном применяется для выращивания кристаллов небольшого диаметра (порядка 2"), с повышенными требованиями по
© Любимова Т.П., Хлыбов О. А., 2003
однородности примеси в кристалле, что необходимо, например, для производства полупроводниковых лазеров.
УВ имеет несколько модификаций (ЬБУВ, УвБ и др.) [1], отличающихся в основном техническими деталями конструкции печи: так, например, в "классическом" УВ ампула движется внутри печи с постоянным распределением температуры, а в процессе УвР движущихся элементов нет, а управление ростом происходит изменением температурного режима нагревателей. Протекающие же в них физические процессы в основном сходны, что позволяет распространять полученные результаты на все семейство.
Магнитное поле является эффективным средством воздействия на поведение проводящей жидкости, позволяющим существенно влиять на структуру и интенсивность течения и, как следствие, на форму фронта кристаллизации и распределение примеси.
Вращающееся магнитное поле (ВМП) [2] позволяет добиваться заметных результатов даже при относительно небольших значениях, порядка 1 мТл, напряженности магнитного поля, что существенно для промышленного производства. ВМП увеличивает общую интенсивность течения в расплаве, а это улучшает гомогенизацию примеси в расплаве и уменьшает прогиб фронта кристаллизации.
1. ПОСТАНОВКА ЗАДАЧИ. ОПРЕДЕЛЯЮЩИЕ УРАВНЕНИЯ И ГРАНИЧНЫЕ УСЛОВИЯ
Модель печи для выращивания полупроводниковых кристаллов по вертикальному методу Бриджмена представляет собой цилиндр с твердыми границами. Толщина границ принята равными нулю. Цилиндр заполнен двумя фазами одного вещества - твердой внизу и жидкой вверху; вещество содержит пассивную примесь. Поверхность раздела фаз в общем случае искривлена (рис. 1).
Формирование теплового режима внутри ампулы и управление процессом роста осуществляется заданием зависящего от времени температурного профиля на боковой поверхности цилиндра. Температура на боковой поверхности задается исходя из линейного закона теплоотдачи. Температура на верхней и нижней границах цилиндра зафиксирована и соответствует в каждый момент времени значению внешнего распределения температуры на соответствующей высоте. На фронте кристаллизации ставится условие Стефана, учитывающее разные теплопроводности фаз, теплоту фазового перехода и движение фронта.
Рис. 1. Вертикальный метод Бриджмена. Модель печи
Система уравнений конвекции в двухкомпонентной среде в приближении Буссинеска при наличии ВМП имеет вид:
ЭК / - \ - 1 - -
— + (УУ)У = -~ГУР +УА¥ - ёРТт + Е1 , (1.1)
УУ = 0, (1.2)
+ {^У)Тт = ХшАГт , (1.3)
1Н =Хс АТ, (1.4)
Эс
— +(УУ)С = ВАС . (1.5)
Здесь индексы с и т означают принадлежность к кристаллу и расплаву, соответственно, Ё1 - сила Лоренца; остальные обозначения - обычные.
Внешнее магнитное поле предполагается однородным, аксиально-симметричным и вращающимся с высокой частотой, что дает возможность ввести скалярный магнитный потенциал / и произвести осреднение силы Лоренца:
р1 = I1 -1 !][й х-], (1-6)
А/-4 = 0. (1.7)
г
При выводе (1.6) и (1.7) влияние индуцированного течением магнитного поля не учитывается в силу предположения о слабости естественной конвекции. Число полюсов магнита принято равным единице [3].
Граничные условия для температуры
На боковой поверхности считается выполненным линейный закон теплоотдачи:
Эп
к
а
ЭТ
Г к
(Т Г- Т ¥ (^ о). (1.8)
Здесь ка - коэффициент теплоотдачи материала ампулы, к1 -коэффициент теплопроводности г-й фазы, Т¥ (г, /) - заданное внешнее распределение температуры, в общем случае зависящее от времени.
На верхней и нижней границах поддерживаются фиксированные температуры, соответствующие внешним температурам на соответствующих высотах:
- для верхней границы: Т|Г - Т¥ (Н,/);
- для нижней границы: Т|г - Т¥ (0, ^ .
На фронте кристаллизации ставится условие Стефана:
ЭТ Кс ЭИ
ЭТ -к„ — Эп
ш
с
+ - 0 СРАТ х
(1.9)
учитывающее движение фронта (У{- скорость перемещения фронта кристаллизации), разные теплопроводности фаз и теплоту фазового перехода АН .
Граничные условия для магнитного поля
Для скалярного потенциала магнитного поля считаются выполненными следующие граничные условия:
- на оси: /\ — 0 ;
■> 1г—0 ’
т
- 0 ;
Э/
- на боковой поверхности: —
Эп
на верхней и нижней границах:
Эп
При выводе граничных условий использовано условие (] ■ п )| - 0, означающее отсутствие нормальных электрических токов на границе (ампула из неэлектропроводящего материала).
Граничные условия для скорости Условие прилипания на твердых границах: к|г - 0 .
Здесь скорость движения фронта кристаллизации принимается много меньшей характерной скорости движения жидкости и потому не учитывается.
ЭК
На оси:
Эп
- 0.
Граничные условия для концентрации примеси
Для концентрации примеси считаются выполненными следующие граничные условия:
ЭС
- на оси, верхней и боковой границах: —
Эп
- на фронте кристаллизации:
ЭС Эп
-0;
+ С(1 - кV - 0
(110)
Здесь к - коэффициент сегрегации, отвечающий за захват примеси кристаллом ( к -я часть примеси по массе остается в кристалле, тогда как (1 - к) - выбрасывается обратно в расплав).
ассмотрение ограничивается случаем осесимметричного решения. В этом случае удобно ввести завихренность р и функцию тока
У:
Р-
Эх Эг
V -1 У
г г Ъх
(1.11)
(1.12)
Г
- г
Г
Г
Г
к =-1 у
2 г дг
(1.13)
Система уравнений (1.1) - (1.7) в терминах функции тока, завихренности и азимутальной скорости Уф , для которой в дальнейшем используется обозначение w , принимает вид:
д Ф _ 1 ( дфду дфду +фду'Л + w дw
д/ г ^ д2 дг дг д2 г д2 ) г д2
I . фФ ёЬ дТ
+,\ ф_г.)- г_
д2у д2у 1 ду _
гф---------------------------------------------------------------^-т- +-—_ 0 (1.15)
^ дг2 дг 2 г дг
дw 1 (дyдw дyдw ду^ ( , w
— _-1 —-------- ----| + п| Aw- —
д/ г ^ дг д2 д2 дг д2 ) + оОВ2 ( _д/
д/ г ^ дг д2 д2 дг )
(1.14)
(1.16)
дТт 1 (дудТт дудТт .
^7 _~\ -Т- ^ _^ I + Ст АТт (1.17)
д/ г ^ дг д2 д2 дг
Щ- _С АТС (1.18)
дС 1 (ду дС ду дС ^
-_-\ ------| + £АС (1.19)
Af-^ _ 0 (1.20)
г
2. МЕТОД РЕШЕНИЯ
При решении задач численного моделирования процессов выращивания кристаллов часто применяется так называемый квазиста-тический подход, при котором фронт кристаллизации предполагается неподвижным, а эффекты, связанные с движением границ, моделируются с помощью граничных условий специального вида. В этом случае задачу можно считать стационарной, что часто позволяет упростить и ускорить ее решение. Однако реальные процессы,
происходящие в расплаве, часто бывают нестационарными. К примеру, существуют процессы, в которых прогиб фронта кристаллизации никогда не выходит на стационарное значение, а имеет тенденцию к увеличению до окончания процесса. В методах выращивания кристаллов, в которых имеются свободные поверхности (например метод жидкой зоны), конвективное течение часто теряет устойчивость и становится нестационарным при достаточно малых значениях определяющих параметров. Поэтому использование ква-зистатического подхода серьезно ограничивает доступный для исследования диапазон параметров. В данной работе используется другой, нестационарный подход, который позволяет обойти эти ограничения.
Система уравнений (1.14) - (1.20) решалась численно методом конечных разностей, по явной схеме с использованием центральных разностей. В физической системе координат (ФСК) сетка криволинейная, не ортогональная. Для спрямления фронта кристаллизации и получения прямоугольной ортогональной сетки в вычислительной системе координат (ВСК) применялось следующее преобразование координат:
Здесь величины со штрихом относятся к ФСК, без штриха - к ВСК; с и т обозначают принадлежность к области кристалла и расплава, соответственно, Е'(г) - функция, определяющая форму фронта в ФСК, Е - 7-координата спрямленного фронта в ВСК, Н
- высота ампулы. Аналогичным образом преобразуются производные в уравнениях и граничных условиях.
Расчет нестационарного процесса роста кристалла производился в несколько этапов:
• Решение чисто теплопроводной задачи во всей ампуле в отсутствие движения нагревателя для определения теплопроводного распределения температуры и формы фронта.
• Решение полной задачи с учетом конвекции в отсутствие движения нагревателя для определения полей температуры, мери-
(2.1)
(2.2)
(2.3)
диональных и азимутальной компонент скорости и формы фронта при фиксированном внешнем распределении температуры Т¥.
• Моделирование процесса выращивания кристалла при нестационарном нагреве с учетом движения нагревателя.
Алгоритм вычислений на каждом шаге по времени включал:
• Вычисление новых полей завихренности, азимутальной скорости, температуры и концентрации из (1.14) и (1.16) - (1.19).
• Нахождение полей функции тока (1.15) и магнитного потенциала (1.20) итерационным методом.
• Определение нового положения и формы фронта кристаллизации Р '(г) по рассчитанному распределению температуры.
• Пересчет преобразования координат на основе полученного Р'(г) .
Таблица 1
Высота Н 5 см
Радиус Я 0.75 см
Начальное положение фронта ¿0 2.5 см
Скорость роста V 0 4 •10-4 с^ с
Температура плавления Тт 937.4 0С
Вертикальный перепад температуры АТ 150 0С
Плотность (кристалл, расплав) Р 5.5 г смъ
Теплопроводность (кристалл) кс 0.17 Вт/см • К
Теплопроводность (расплав) к 0.39 Вт/ см • К
Теплоемкость (кристалл, расплав) СР 0.39 Дж/г • К
Удельная теплота фазового перехода АН 460 Дж/г
Коэффициент теплового расширения Ь 5 •10-4 к-1
Кинематическая вязкость V 1.3 -10-3 см1 / с
Удельная электропроводность о 2.37-106 (Ом•м)1
Частота вращения поля О 50 с-1
Конвективные течения. Сборник научных трудов
3. ЧИСЛЕННЫЕ РЕЗУЛЬТАТЫ
Численные эксперименты проводились для процесса выращивания полупроводникового кристалла ваіве (германий с примесью галлия). Использовалось линейное распределение температуры Т¥, которое рассчитывалось исходя из общей высоты ампулы Н , перепада температуры АТ и зависящего от времени положения нагревателя 2 (/). Основные параметры приведены в табл. 1.
Влияние ВМП на форму фронта кристаллизации
Рис. 2 иллюстрирует влияние ВМП на прогиб фронта кристаллизации. Представлена эволюция прогиба фронта в процессе роста кристалла для трех различных интенсивностей магнитного поля. В момент времени X происходит запуск конвекции (до этого момента учитывается только диффузионный механизм теплопереноса); момент У соответствует установившемуся режиму при неподвижном нагревателе; момент 2 - выходу на плато - установлению квазиста-ционарного режима роста.
сек
ХУ г
Рис. 2. Зависимость прогиба фронта кристаллизации от времени
При всех значениях параметров, для которых проводились расчеты, фронт кристаллизации получается вогнутым. Искривление фронта происходит даже в отсутствие конвективного течения и движения нагревателя; причиной этому служит различие теплопроводностей фаз и теплота фазового перехода [4].
Прогиб фронта максимален для чисто диффузионного режима теплопроводности и уменьшается с появлением конвекции. Движение нагревателя приводит к некоторому увеличению прогиба. ВМП приводит к уменьшению кривизны фронта. При достижении некоторой критической величины магнитного поля (в данном случае ~0.7 мТл) азимутальное течение теряет устойчивость, в результате чего положение и форма фронта начинают изменяться хаотически.
Влияние ВМП на структуру и интенсивность течения
На рис. 3 приведены изолинии функции тока в меридиональном сечении в момент времени 2 (рис. 2) для трех различных значений интенсивности поля. В отсутствие магнитного поля течение имеет одновихревую структуру; центр вихря располагается вблизи фронта; вихрь вращается по часовой стрелке (рис. 3, слева; В = 0, у =-0.014; здесь и далее показана только правая часть сечения ампулы).
Рис. 3. Изолинии функции тока ¥ при различных интенсивностях магнитного поля: В = 0 , 0.7 и 0.8 мТл
Магнитное поле приводит к образованию в верхней части ампулы второго вихря, вращающегося против часовой стрелки (рис. 3, в центре; В = 0.7 мТл, у = -0.015 + 0.024). При увеличении магнитного поля верхний вихрь распадается на три более мелких (рис. 3, справа; В = 0.8 мТл, у = -0.015 + 0.030).
ш г 5.00- 4.50- (И г 5.00- 4.50- (Ш-
'Ж
- 3.50- ЛЩУ - 3.50- Іі'Ш-
•- 1 - 3.00- - 3.00- /Щ
- 2.50- і і - 2.50- Шт-
- 2.00- - 2.00-
- 1.00- - 1.00-
- 0.50- - 0.50-
- 0.00- - 0.00-
0.00 0.50 0.00 0.50
Рис. 4. Эволюция функции тока при В = 1.0 мТл; А/ = 10 с
При дальнейшем усилении магнитного поля течение становится нестационарным; появляются характерные вихри Тейлора [5].
На рис. 4 показана перестройка структуры течения в меридиональном сечении со временем при В = 1.0 мТл. Период между кадрами 10 сек. На рис. 5 представлены соответствующие азимутальные скорости с максимальным значением ~ 0.5 см/с. Течения показаны на момент выхода на квазистационарный режим роста (точка 2 на рис. 2).
Влияние ВМП на распределение примеси
В силу больших чисел Шмидта для расплава ва:ве, конвективный механизм переноса примеси преобладает над диффузионным - наблюдается "вмороженность" изолиний концентрации примеси в расплаве в линии тока (рис. 6).
0.00 0.50
0.00 0.50
0.00 0.50
Рис. 6. Распределение примеси при различных интенсивностях магнитного поля: В = 0 , 0.7 и 0.8 мТл
ВМП увеличивает общую интенсивность меридионального течения, что, в свою очередь, усиливает выброс примеси в расплав - что можно видеть по увеличению числа изолиний концентрации на рис. 6. (Изолинии концентрации приведены для момента
установления квазистационарного режима роста - точка 2 на рис. 2.)
5.00
5.00
5.00
4.50
4.50
4.50
4.00
4.00
4.00
3.50
3.50
3.50
3.00
3.00
3.00
2.50
2.50
2.50
2.00
2.00
2.00
1.50
1.50
1.50
1.00
1.00
1.00
□ .50
□.50
0.50
0.00
0.00
0.00
Конвективные течения. Сборник научных трудов БИБЛИОГРАФИЧЕСКИЕ ССЫЛКИ
1.Rudolph P., Matsumoto F., Fukuda T. Studies on interface curvature during vertical Bridgman growth of InP in a flat-bottom container // J. of Crys. Growth. 1996. V. 158. P. 43-48.
2. Dold P., Benz K.W. Rotating magnetic fields: fluid flow and crystal growth applications // Prog. In Crys. Growth and Char. of Mat. 1999. P. 7-38.
3.Горбунов Л.А., Колевзон В.Л. Течение проводящей жидкости во вращающемся магнитном поле // Магнитная гидродинамика. 1992. № 4. С. 69-74.
4.Lan C.W., Ting C.C. Numerical investigation on the batch characteristics of liquid encapsulated vertical Bridgman crystal growth // J. of Crys. Growth. 1995. V. 149. P. 175-186.
5. On the stability of rotating MHD flows / Marty P., Witkowski L.M., Trombetta P., etc. // Transfer Phenomena in Magnetohydrodynamic and Electroconducting Flows. Vol. 51. Fluid Mechanics and its Applications / Eds.: A. Alemany, P. Marty and J.P. Thibault. Kluwer Academic Publishers, 1998. P. 327-343.
NUMERICAL SIMULATION OF GA:GE VERTICAL BRIDGMAN CRYSTAL GROWTH UNDER THE INFLUENCE OF ROTATING MAGNETIC FIELD
T.P. Lubimova, O.A. Khlybov
Abstract. The results of 2D transient analysis of Ga:Ge vertical Bridgman crystal growth in flat-bottom ampoule are presented. Solid-liquid interface shape temporal evolution is considered.
The effect of rotating magnetic field on interface shape, flow patterns and intensity of both meridional and azimuthal flows, and dopant segregation is studied. Transient flow patterns for the case when azimuthal Taylor - Couette flow becomes unstable due to over-critical magnetic field intensity, along with the development of Taylor vortices, are shown.
