Численное моделирование процессов тепло- и массообмена при выращивании монокристаллов методом Чохральского Текст научной статьи по специальности «Физика»

УДК 532.534

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ПРОЦЕССОВ ТЕПЛО- И МАССООБМЕНА ПРИ ВЫРАЩИВАНИИ МОНОКРИСТАЛЛОВ МЕТОДОМ ЧОХРАЛЬСКОГО

ТЕВЯШЕВ А.Д., СУЗДАЛЬ В. С, БОРОДАВКО Ю.М, ПЕЛИПЕЦА.А.

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

Кристаллы широко применяются в элементной базе ЭВМ, лазерной технике, средствах связи. Определяющую роль в формировании монокристаллов играют гидромеханические процессы и процессы тепло- и массопереноса. В [3] дано обоснование математической модели процесса выращивания монокристаллов методом Чохральского. Реализация этой модели на ЭВМ приведена ниже.

1. Схема установки

Рассмотрим установку для выращивании монокристаллов (рис. 1).

Рис.1. Схема установки для выращивания кристаллов методом Чохральского: Rk — радиус кристалла; Rj — радиус тигля; H — глубина расплава; Q к — угловая скорость вращения кристалла; Q.j — угловая скорость вращения тигля

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

2. Математическая модель

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

В этой работе используется формулировка уравнений Навье-Стокса в переменных ^ , а> , W при численном решении методом конечных разностей (МКР) в односвязной области простой геометрической формы с прямолинейными границами:

дю ттдю т,дю Urn 1 dm2)

--+ U---+ V------------*--<■ =

dt dr dz r r dz

1 Г Gr de Grd dC (1)

Re ^ r2 J Re2 dr Re2 dr

d 2y 1 dy + d 2y Or2 r dr Qz2

(2)

dW TTdW TdW UW

---+ U----+ V----+----

dt dr dz r

ReW

2,

(3)

dU

где Ю = —

dz

dV

dr

U =1 ^

r dz

V = -1 ^ . r dr

Здесь U , V , W — компоненты вектора скорости;

(О — вихрь; ц/ — функция тока; Д

52 1 5 52

—2 ^----1---2

dr2 r dr dz2

—лапласиан скалярной функции, в качестве характерных масштабов выбраны скорость вращения

Qк и радиус Rk кристалла.

Совместно с системой (1)-(3) решаются уравнения переноса тепла (4) и примеси (5):

dd TTde T,dd 1

— + U— + V— Ав

dt dr dz RePr

dC dC dC 1

+ U + V = -AC

dt dr dz RePrd

(4)

(5)

Начальные условия: ^ = 0 ; W = 0; в = 0; C = 0 .

12

РИ, 2002, № 1

Граничные условия:

G, \w = 0,= 0, W = — r,U = 0,V = 0,в = виC = 1; 1 dz QK 1

G2 \w = 0,^ = 0,W = -^T-RT,U = 0,V = 0,0=02,C = 1; 2 dr Q.KT 2

^ „ Mn d0 Mnd dC dW

G3 :w = 0,a =----+-----d--,----= 0,

RePr dr Re Prd dr dz

U =Mnn_ M+Mn^ 5C ,V = 0,0 =e ^ = 0;

RePr dr Re Prd dr dz

дш dC

G4 :p = 0,-^- = 0,W = r,U = 0,V = 0, 0 = 04,— = BidC; dz dz

r)V г)й r)C

G5 :w = 0,® = 0,W = 0,U = 0,— = 0,— = 0,— = 0 ,

dr dr dr

где G, = {(r,z)|0 < r < Rt,z = 0} — дно тигля; G2 = {(r , z)|r = Rt ,0 < z < И}—боковая стенка тигля; G3 = {(r, z~)RK < r < Rt , z = и}—поверхность расплава; G4 = {(r, z)|0 < r < RK, z = и} — фронт кристаллизации; G5 = {(r,z'y = 0,0 < z < и} — ось симметрии.

T - TK

Параметры системы: 0 =

AT

безразмерная

величина температуры; AT = Tt - Tk — максимальная разница температур в расплаве; Tk , Tt — максимальная и минимальная величина температур на поверхности кристалла и стенках тигля соответственно; T) — температура окружающей среды;

.0 = 0— безразмерная температура окружаю-

щей среды; Ck , Ct — концентрация примеси в кристалле и на стенках соответственно; AC = Ck - Ct — максимальная разность концент-

рации примеси; Re =

Q KR

к^к_

v

— число Рейнольдса;

з

Gr =-------— AT — число Грасгофа (тепловое);

2

Gr =

gPdRK

2

AC — число Грасгофа (концентраци-

vv онное); Pr =--число Прандтля; Prd = Sc = ^ —

число Прандтля (концентрационное);

Ра oRK

Mn = —т----AT — число Марангони (тепловое);

pva

Mn P^daRK

Mnd =

pvD

AC — число Марангони (концент-

aR

рационное); Bi = —к — число Биота (теговое);

a

(l - k )ufRK

Bid =---d------_ число Биота (концентрацион-

РИ, 2002, № 1

ное); Ped = RePrd — число Пелле; Ra = Gr Pr — число Рэлея (тепловое); Rad = Gr Prd — число Рэлея (концентрационное); g — ускорение силы тяжести; v, d , a — коэффициенты кинематической вязкости, диффузии и температуропроводности; Ud — фиксированная скорость перемещения фронта кристаллизации; Р- , Pd — коэффициенты объемного расширения — теплового и концентрационного; Рат , Pad — коэффициенты поверхностного расширения — теплового и концентрационного; а — коэффициент теплообмена на поверхности расплава; k — коэффициент распределения примеси на фронте кристаллизации.

3. Метод решения

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

Для аппроксимации дифференциальных уравнений введем пространственно-временную сетку с координатами ri = ih , zt = jl, ті = nr , где h , l — шаги сетки по координатам r , z соответственно; г - шаг по времени, i = 0,M , j = 0, N, n = 0, K .

Введем следующее обозначение: p(ih, jl, nr)= pfj . Производные по пространственным переменным будем аппроксимировать конечными разностями:

dp Фі+1, j _ Pi-1, j d 2p ~ Pi+l,j 2Pi,j + Pi-l, j

dr 2h ’dr2 h2

dp Pi, j+1 -Pi, j-1 d 2p Pi, j +1 - 2Pi, j +Pi, j-1

~dz 2l ’ dz2 l2

Производную по времени заменим разностным отношением “вперед” в виде

л n+1 n

dp Pi, j Pi, j

dz z

Условие устойчивости явной схемы z <

ch2

где c

— параметр, зависящий от числа Рейнольдса и убывающий от 1 до 0,1 при увеличении этого числа от 0 до 600.

4. Анализ результатов

С помощью построенной математической модели проведем исследование поведения линий тока в зависимости от скорости вращения тигля и кристалла, размеров выращиваемого кристалла по отношению к размерам тигля, от соотношения высоты тигля и его ширины, числа Рейнольдса Re .

Параметры системы:

R— — 1; Rk ~ 0,5; ^к — 1; ^т — —0,5; Re = 300; Gr = 1-105; GrD = 1 • 103; Pr = 0,02 ; PrD = 0,3; Mn = 5,1 • 104 ;

13

Bi = 0,5 ; BiD =-15,6, m=21, n=21.

Граничные условия: T1 = T2 = 1, T3 = 0,5, T4 = 0 .

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

Рис.2. Линии уровня функции тока: Q к = 1, О t =-0,5

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

Рис. 3. Линии уровня температуры: Q к = 1,0 т = -0,5

20

15

10

5 •

5 10 15 20

Рис.4. Линии уровня концентрации: Q к = 1, ^т = -0,5

При вращении только кристалла или при отсутствии вращения наблюдается одновихревое течение (рис. 5, 6).

Рис. 5. Линии уровня функции тока: = 0 , D-т = 0

Рис.6. Линии уровня функции тока: Q к = 1, ^т = 0

14

РИ, 2002, № 1

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

20

15

10

5

0 •

5 7.5 10 12.5 15 17.5 20

Рис. 7. Линии уровня функции тока:

Qк = 1,0т = -0,5 , Re = 688

Рис. 8. Линии уровня функции тока: Пк = 1, Q.T = -0,5 , Re = 1,57 -104

При изменении радиуса кристалла (рис. 9,10) изменяются размеры и интенсивность подкристальной циркуляции: чем больше радиус кристалла, тем крупнее вихрь и выше скорость движения расплава. Следует отметить, что при увеличении подкристальной циркуляции возрастает угроза разъединения кристалла и расплава.

Рис. 9. Линии уровня функции тока: О к = 1,0 т = —0,5 , Rk = 0,2

Рис. 10. Линии уровня функции тока:

О к = 1,0 т = —0,5 , Rk = 0,7

На рис. 11, 12 представлены линии уровня функции тока для различного числа точек разбиения области. При увеличении числа точек разбиения повышается точность решения, и линии уровня становятся более гладкими (рис. 12).

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

(Ит / Rt) и 1, О к ~ 1, ^ т ~ “0,5 , Rk и 0,5 ; подвод тепла можно осуществлять через нижнюю и боковую стенки.

20

15

10

РИ, 2002, № 1

15

5. Заключение

20 •

15

10 •

5

0 •

0 5 10 15 20

Рис. 11. Линии уровня функции тока:

QK = 1, Qf = “0,5 , m = 21, n = 21

30 •

25 • •

20 • . '

15 •

10 •

5 •

0 •

0 5 10 15 20 25 30

Рис. 12. Линии уровня функции тока:

QK = 1, ^т = —0,5 , m = 30 , n = 30

УДК 621.382:537.5

ВОЗМОЖНОСТИ ЭКСПЛУАТАЦИИ ЧИСЛЕННОЙ ДИНАМИЧЕСКОЙ МОДЕЛИ ЭЛЕКТРОННО-ЛУЧЕВЫХ ПРИБОРОВ

ГЛУМОВА М.В., ВОРОБЬЕВ М.Д.___________

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

Научно-исследовательское проектирование электронных приборов требует наличия надежных программных кодов. В качестве одного из таких программных продуктов авторы предлагают использовать разработанную численную динамическую модель электронно-лучевых приборов (ЭЛП).

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

Литература: 1. Математическое моделирование конвективного тепломассообмена на основе уравнений Навье-Стокса / В.И. Полежаев, А.В. Бунэ, Н. А Верезуб и др. М.: Наука, 1987. 272 с. 2. Полежаев В.И., Пасконов В.М., Чудов Л.А. Численное моделирование процессов тепло- и массообмена. М.: Наука, 1984. 288с. 3. Математическое моделирование физических процессов при выращивании монокристаллов методом Чохральского / А.Д. Тевяшев, В.С. Суздаль, Ю.М. Бородавко, А.А Пелипец // Радиоэлектроника и информатика. 2001. № 4. С.33-43. 4. Алексеенко С.В. Математическое моделирование процессов тепло- и массообмена при выращивании кристаллов по методу Чохральского // Вісн. Дніпропетр. ун-ту. 1999. Вип. 2. Т.1. С. 100-104.

Поступила в редколлегию 19.11.2001

Рецензент: д-р физ.-мат. наук, проф. Гордиенко Ю.Е.

Тевяшев Андрей Дмитриевич, д-р техн. наук, профессор, зав. каф. ПМ ХНУРЭ. Научные интересы: системный анализ и теория оптимального стохастического управления. Адрес: Украина, 61166, Харьков, пр. Ленина, 14.

Суздаль Виктор Семенович, канд. техн. наук, зав. отделом НТК “Институт монокристаллов”. Научные интересы: системы управления технологическими процессами получения монокристаллов. Адрес: Украина, 61001, Харьков, пр. Ленина, 60.

Бородавко Юрий Михайлович, канд. техн. наук, доцент кафедры ПМ ХНУРЭ. Научные интересы: математическое моделирование физических процессов. Адрес: Украина, 61166, Харьков, пр. Ленина, 14.

Пелипец Андрей Александрович, аспирант кафедры ПМ ХНУРЭ. Научные интересы: математическое моделирование процессов выращивания монокристаллов из расплава. Адрес: Украина, 61166, Харьков, пр. Ленина, 14.

Целью ее создания явилось моделирование физических процессов в ЭЛП с максимальным приближением к реальным условиям функционирования этих устройств. Разработанная модель — инструмент проведения численных экспериментов для исследования как установившихся, так и переходных процессов в объеме ЭЛП, с возможностью задания начальных и граничных условий при разных режимах работы устройств с учетом полей пространственного заряда для пучков различной интенсивности, со статистическим моделированием процесса эмиссии и учетом прикатодных явлений, с возможностями исследования флуктуацион-ных процессов.

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

16

РИ, 2002, № 1