CC BY
73
УДК 532.534
МАТЕМАТИЧЕСКИЕ МОДЕЛИ ФИЗИЧЕСКИХ ПРОЦЕССОВ ПРИ ВЫРАЩИВАНИИ МОНОКРИСТАЛЛОВ МЕТОДОМ ЧОХРАЛЬСКОГО
ТЕВЯШЕВ А.Д., СУЗДАЛЬ В.С.,
БОРОДАВКО Ю.М, ПЕЛИПЕЦА.А._____________
Рассматриваются уравнения гидромеханики и тепломассопереноса, составляющие единую математическую модель процесса выращивания монокристаллов из расплава методом Чохральского. Обсуждаются приближения этой модели, в той или иной степени адекватные реальным условиям роста монокристаллов.
1. Введение
Определяющую роль в формировании обьемных монокристаллов (МК), используемых в электронике, играют гидромеханические процессы и процессы тепло- и массопереноса [1,2,4,5]. Исследование этих физических процессов до последнего времени проводилось в основном экспериментальными методами. Однако современное приборостроение предьявляет все более высокие и разнообразные требования к свойствам и качеству МК. Это резко усложняет экспериментальные исследования, их стоимость и трудоемкость. В связи с этим возникает необходимость использовать в исследованиях, направленных на совершенствование и оптимизацию способов получения МК, методы математического моделирования.
Теоретические исследования призваны обеспечить достаточно точные и надежные рекомендации, необходимые для проектирования оптимальных конструкций и режимов работы установок по выращиванию МК. Это требует рассмотрения большого числа взаимосвязанных процессов, учета де -талей геометрии ростовых установок, достаточно точного описания свойств изучаемых сред. Поэтому при построении математических моделей необходимо, чтобы они, с одной стороны, были достаточно полны и могли обеспечить требуемую точность и надежность результатов, а с другой — реализуемы на существующих ЭВМ.
2. Метод Чохральского
Выращивание МК по методу Чохральского состоит в том, что монокристалл медленно вытягивается на затравку со свободной поверхности расплава и одновременно вращается вокруг своей вертикальной оси (рисунок): 1 — затравка, 2 — кристалл, 3 — фронт кристаллизации, 4 — расплав, 5 — тигель, R — радиус кристалла, h — высота фронта кристаллизации, H — уровень расплава в тигле, V — скорость вытягивания кристалла, q — угловая скорость вращения кристалла, zOr — система координат.
В результате вращения получают симметричные кристаллы совершенной структуры. Одновременно с выращиванием монокристаллов происходит эф-
фективная очистка вещества от примесей. Это обусловило широкое применение данного метода для выращивания МК высокой степени чистоты
Расплав, находящийся в цилиндрическом или коническом тигле, подогревается обычным или индукционным нагревателем, причем температура расплава поддерживается несколько выше точки плавления данного вещества. Монокристаллическая затравка определенной ориентации крепится к вытягивающему механизму (кристалл од ержателю) и опускается в расплав сверху для затравливания, а затем, после некоторого промежутка времени, начинает медленно извлекаться при вращении вокруг своей оси. При подьеме монокристалл, смачиваясь, увлекает за собой часть жидкости, кристаллизующейся по мере подьема кристаллодержателя. Скорость извлечения последнего из расплава равна скорости роста кристалла, которая в значительной степени зависит от отвода тепла через затравку.
Необходимо отметить, что расплав можно представить как вязкую неоднородную по составу жидкость. Уравнения гидродинамики, тепло - и массопереноса в такой жидкости можно вывести на основе уравнений гидродинамики идеальной жидкости.
Схема выращивания кристаллов методом Чохральского
3. Уравнения гидродинамики идеальной жидкости
Математическое описание состояния движущейся жидкости осуществляется с помощью функций, определяющих распределение скорости жидкости:
и = ії(х, у, z, t)
и каких-либо ее двух термодинамических величин, например давления px, у, z,t) и плотности
р{ X У, z, 0.
Как известно [3], все термодинамические величины определяются по значениям каких-либо двух из них с помощью уравнения состояния вещества; поэтому задание пяти величин: трех компонент скорости и , давления p и плотности р полностью определяет состояние движущейся жидкости.
К числу наиболее важных гидродинамических уравнений относится уравнение непрерывности.
РИ, 2001, № 4
33
3.1. Уравнение непрерывности
Это уравнение представляет закон сохранения вещества в гидродинамике.
Рассмотрим некоторый объем V0 пространства.
Количество (масса) жидкости в этом объеме | pdV , где р — плотность жидкости. V
Через элемент df поверхности, ограничивающей объем V0, в единицу времени протекает количество podf жидкости; вектор df по абсолютной величине равен площади элемента поверхности и направлен по нормали к ней.
Полное количество жидкости, вытекающей в единицу времени из объема V0 , есть, следовательно,
jp»df, (1)
где интегрирование производится по всей замкнутой поверхности, охватывающей объем Vo .
С другой стороны, уменьшение количества жидкости в объеме Vo можно записать в виде
Приравнивая (1) и (2), получаем:
Поскольку
(2)
(3)
(3) преобразуем к виду
^\^ + div{pv)jc/P=0
и соответственно к виду
(4)
(5)
div[pu) = 0
(6)
— уравнение непрерывности.
3.2. Уравнение Эйлера
Это уравнение, также как и уравнение непрерывности, относится к числу основных уравнений гидродинамики.
Выделим в жидкости некоторый объем. Полная сила, действующая на выделенный объем жидкости, равна интегралу
Преобразуя его в интеграл по объему, имеем:
-dVgradp. (9)
Иными словами — на единицу объема действует сила
-gradp . (10)
Тогда уравнение движения элемента объема жидкости будет иметь вид
dd
Р— = ~gradp . (11)
Следует заметить, что изменение do скорости элемента объема жидкости в течение времени dt складывается из двух частей:
а) из изменения скорости в данной точке пространства в течение времени dt;
б) из разности скоростей (в один и тот же момент времени) в двух точках, разделенных расстоянием df, пройденным рассматриваемым элементом объема жидкости в течение времени dt .
Первая из этих частей равна
до ,
---dt
dt :
(12)
вторая
dx — + dy — + dz — = (drV)o
dx dy dz .
Тогда (11) преобразуется к виду
do 1 ,
----VI oV p =-gradp
dt p
— уравнение Эйлера — уравнение движения жидкости.
(13)
(14)
Если жидкость находится в поле силы тяжести, то на каждую единицу ее объема действует еще сила pg . Тогда (14) приобретает вид
do 1 ^ _
— + (oV)o =-----Vp + g
dt p
(15)
При выводе уравнений движения жидкости не учитывались процессы диссипации энергии, имеющие место в текущей жидкости вследствие внутреннего трения (вязкости) в ней и теплообмена между различными ее участками. О таком движении жидкости говорят как о движении идеальной жидкости.
Отсутствие теплообмена между отдельными участками жидкости (а также между жидкостью и соприкасающимися с нею окружающими телами) означает, что движение происходит адиабатически. Движение идеальной жидкости—адиабатическое движение.
При адиабатическом движении энтропия каждого участка жидкости остается постоянной при перемещении последней в пространстве. Обозначая посредством 5 энтропию, отнесенную к единице массы жидкости, имеем
ds
dt
= 0
для адиабатичности движения, или
----Ъ ograds = 0
dt
(16)
(17)
34
РИ, 2001, № 4
Это общее уравнение, выражающее адиабатичность движения идеальной жидкости. Используя (6), представим (17) в виде
д
— (ps) + div(psu) = 0 , (18)
где psu — плотность потока энтропии.
Обычно уравнение адиабатичности принимает гораздо более простую форму. Если, как это обычно имеет место, в некоторый начальный момент времени энтропия одинакова во всех точках объема жидкости, то она останется везде одинаковой и неизменной со временем и при дальнейшем движении жидкости. Тогда
s = const. (19)
Такое движение, при котором выполняется условие (19), называется изэнтропическим.
3.3. Уравнение Эйлера и уравнение непрерывности для несжимаемой жидкости
В ряде случаев при течении расплава (жидкости) плотность можно считать не изменяющейся, т.е. постоянной вдоль всего объема жидкости в течение всего времени движения. Другими словами, в этих случаях при движении не происходит заметных сжатий или расширений жидкости. О таком движении говорят как о движении несжимаемой жидкости.
Общие уравнения гидродинамики сильно упрощаются при применении их к несжимаемой жидкости. Уравнение Эйлера не меняет своего вида, если положить в нем р = const, за исключением только того, что в уравнении (15) можно внести р под знак градиента:
f- + №У> = -v£+if. (20)
Зато уравнение непрерывности принимает при р = const простой вид:
divu = 0 (21)
—уравнение непрерывности для несжимаемой жидкости.
4. Вязкая жидкость: уравнение Навье-Стокса
Для того чтобы получить уравнение Навье-Стокса, описывающее движение вязкой жидкости, необходимо ввести дополнительные члены в уравнение движения идеальной жидкости. Что касается уравнения непрерывности, то, как следует из самого его вывода, оно относится в равной мере и к движению вязкой жидкости. Уравнение же Эйлера должно быть изменено.
Будем производить вычисления в тензорных обозначениях. Имеем:
д
¥ ^ =р
duL +д^ dt dt
(23)
Воспользуемся уравнением непрерывности, написав его в виде
t ~k ^ ■
и уравнением Эйлера в форме
Тогда получим
(24)
(25)
д dui dp д
—рог =~РЩ— ог —
dt dxfo cXi dxfo
(p°k)
dp d
H—TTpuiuk-
cXi dxk
(26)
Первый член справа напишем в виде
dp _ s dp
dxi lk dxk
(27)
и найдем окончательно
d d ^
Jt pu-’—t n “ • где тензор Пik определяется как
П ik = psik + P»i»k
(28)
(29)
и называется тензором плотности потока импульса.
Поток импульса, определяемый формулой (29), представляет собой чисто обратимый перенос импульса, связанный с механическим передвижением различных участков жидкости из одного места в другое и с действующими в жидкости силами давления. Вязкость (внутреннее трение) жидкости проявляется в наличии дополнительного переноса импульса из мест с большей в места с меньшей скоростью.
Поэтому уравнение движения вязкой жидкости можно получить, прибавив к “идеальному” потоку импульса (29) дополнительный член crk , определяющий необратимый, “вязкий” перенос импульса в жидкости. Таким образом, тензор плотности потока импульса в вязкой жидкости имеет вид
Пik = psik + PDlDk~ aik = ~aik + pu^k . (30)
Тензор
°ik = - psik +°'ik (31)
назовем тензором напряжений, а a'ik — вязким тензором напряжений.
Уравнение Эйлера можно записать в виде [3]:
d d ^
РЩ = -^пik .
dt
dxk
(22)
где П ik — тензор плотности потока импульса, определяемый следующим образом.
Импульс единицы объема жидкости есть ри .
д ^
Определим скорость его изменения — ри .
Установить общий вид тензора а'ік можно исходя из следующих соображений. Если градиенты скорости не очень велики, то можно считать, что обусловленный вязкостью перенос импульса зависит только от первых производных скорости. Саму
doi
зависимость a’ik от производных можно в том
же приближении считать линейной.
РИ, 2001, № 4
35
Относительно просто можно доказать, что наиболее общим видом тензора второго ранга a'ik является
ґ
°ik = Л
до, док 2
до,
dxk dxi
— —
3 дх
I
+{Slt Ур. (32)
с независящими от скорости коэффициентами і?и(. Величины 7 и С называются коэффициентами вязкости (причем С часто называют второй вязкостью).
Уравнения движения вязкой жидкости можно теперь получить непосредственно путем прибавления
da'ik _
выражения ~^х к правой части уравнения Эйлера
Р
doi
dt
' + °k
до,■
дх,
к
dp
дх.
(33)
Таким образом, получаем
+
(34)
Это есть наиболее общий вид уравнений движения вязкой жидкости. Величины р и ^ являются, вообще говоря, функциями давления и температуры. В общем случае p, т , а потому и р, £ не постоянны вдоль всей жидкости, так что р и С, не могут быть вынесены из-под знака производной.
В большинстве случаев, однако, изменение р и Q вдоль жидкости незначительно, поэтому можно считать их постоянными. Тогда уравнение (34) можно представить в векторном виде:
Р
до
dt
gradp + рАо +|V +~jgraddivo .(35)
Это так называемое уравнение Навье-Стокса (НС). Оно существенно упрощается, если жидкость можно считать несжимаемой. Тогда divo = 0 и последний член справа в (35) исчезает.
Уравнение Н-С для несжимаемой вязкой жидкости будет иметь вид:
до 1 , n . _
—+ (pV)o= — gradp + — Ao. (36)
dt P P '
Далее уравнение (36) и уравнение непрерывности целесообразно представить в цилиндрической системе координат г,в, z и прямоугольной системе координат х, y, z .
Три компоненты уравнения Н-С в цилиндрической системе координат принимают вид:
dor
dt
+(oV)or
49
r
1 3p+,j
р dr 1
r
2 doQ
Г2 ~дв
,(37)
дщ_
dt
r
11 dp ( oq
---+ H koe----^ +
pr do
2 dor і
r2 ~д0/
(38)
doz
dt
+ {oV)oz
1 dp
p dz
+ vAoz
(39)
Л t
где v = — — кинематическая вязкость (p — динамическая вязкость).
Заметим, что операторы (6V) и д определяются формулами
Л 1 д( д А=----1 r— 1 +
5 1 д д
У-— dr r дв + о7 z dz , (40)
д2 д2 д2 1 5 1 д2 д2
r dr У dr) r2 дв2 dz2 dr2 r dr r2 дв2 dz2'
(41)
Существенно отметить, что при преобразовании слагаемых (бУ)б и дб в (36) использованы следующие тождества векторного анализа:
о2
(бУб) = grad+ [rot б, б] ,
Ао = graddivo-rotroto .
Уравнение непрерывности:
1 д / \ 1 dog doz
11
--\ror I+—
r dr r дв
dz
= 0
или
dor or 1 doQ doz
- H-----h -
dr r r дв
dz
= 0
(42)
(43)
(44)
(45)
Соответственно в прямоугольной системе координат три компоненты уравнения Н-С имеют вид:
бо х
- + о х
бо х
dt дх
1 dp
■ + оу
(д2 о у
р дх
- + V
бод
ду д 2 о
■ + о,
dz
бо у
dt
■ + о х
дх
У
2
ду
5о
X +б 2 О х Л
dz
2
(46)
дх
+ о у
у
1 dp Р ^у
в2
+ V
гу
5 2 о
+ о,
бО у
dz
dx2
У +б 2 о у ^
5о,
■ + о х
5о,
dt дх
_ 1 dp
■ + о у
р dz
■ + v
^ 52 о,
бу
5о z
~ду
д2 о
- + о,
dz
бо
2
(47)
dz
дх
2
Z + 3 2 о z Л
ду
2
dz
2
(48)
Уравнение непрерывности:
дох дОу QUz
■ = 0
дх ду dz ' (49)
5. Уравнение переноса тепла в вязкой жидкости
Вывод этого уравнения основан на законе сохранения энергии в вязкой жидкости. Последний можно получить анализируя закон сохранения энергии в идеальной жидкости, выражаемый уравнением:
д_
dt
( 2 А ( 2 \]
ро V рє 2 V У = —div ро о Ъ w 2 V У_
(50)
где є — внутренняя энергия единицы массы жидкости; w — тепловая функция единицы массы жидкости.
z
36
РИ, 2001, № 4
Слева в (50) стоит скорость изменения энергии единицы объема жидкости, справа—дивергенция плотности потока энергии. В вязкой жидкости закон сохранения энергии тоже имеет место. Однако плотность энергии будет иметь иной вид.
( 2 \
Прежде всего помимо потока Ри
\ /
ного с простым переносом массы жидкости при ее движении, имеется еще поток, связанный с процессами внутреннего трения. Этот второй поток выражается вектором (-бет') с компонентами oia'ik .
Действительно, наличие вязкости приводит к появлению потока импульса a'ik ; перенос же импульса всегда связан с переносом энергии, причем поток энергии получается из потока импульса умножением на скорость. Этим, однако, не исчерпываются все дополнительные члены в потоке энергии.
Если температура жидкости не постоянна вдоль ее объема, то наряду с обоими указанными механизмами переноса энергии будет происходить ее перенос из мест с более высокой в места с более низкой температурой. Он не связан с макроскопическим движением и происходит также и в неподвижной жидкости.
■ + w
,связан-
2
Обозначим через q — вектор плотности потока тепла, переносимого посредством теплопроводности. Он связан некоторым образом с изменением температуры вдоль жидкости. Эту зависимость можно написать сразу в тех случаях, когда градиент температуры в жидкости не слишком велик; тогда можно разложить q в ряд по степеням градиента температуры, ограничиваясь первыми членами разложения:
q = -kvt , (51)
где к — постоянный коэффициент теплопроводности, в общем случае — функция температуры и давления.
Таким образом, полная плотность потока энергии в жидкости при наличии вязкости и теплопроводности равна:
ри
2
- + w
-(и ct')-XVT
(52)
Соответственно этому общий закон сохранения энергии выражается уравнением
д_
dt
2 f 2 Л ~
ро ъ рє 2 V У = -div рб О У w 2 V J -(б a')-XVT
.(53)
Представим его в другом виде, используя уравнение непрерывности, уравнение Навье-Стокса и некоторые термодинамические соотношения. Вычислим производную по времени от энергии единицы объема жидкости:
д_
dt
( 2 Р°
2
| о2 до _ до де др
^ = ТТ+“Р1Т+р¥+‘"а • (54)
v /
др
Подставляя сюда — из уравнения непрерывности
до
и — из уравнения Навье-Стокса, получаем РИ, 2001, № 4
д_
dt
( 2 Л 2 2
ро О _
~Y~+рє =~~2~ divP°-p{°v—
da'ik де ^
+ oi — + р----------edivpo .
dxfr dt
дУр +
(55)
Воспользуемся теперь термодинамическим соотношением
откуда
de = Tds - pdV = Tds H—— dp 2 (56)
P
Зє Dt = T — +44 = T *-4- divl pH) dt p'2 dt dt p2 ' (57)
Тогда, вводя тепловую функцию w = в н—, нахо-
Р
дим
д_
dt
( 2 ри_
2
X (
+ рє
.2 )
w +-
2
2
div( po)-p(oV) — -
_ ds da'ik
-о Vp + pT — + »i-A
Dt OXfc
(58)
Далее из термодинамического соотношения
dp
dw = Tds +— имеем:
P
Vp = pVw - pTVs . (59)
Последний член в правой стороне равенства (58) можно записать в виде
д / , \ , дОі , дОі
4~р— = ^\Рі°ік)= divlucr .(60)
UXfc иХк vXk vXk ^ ^
Подставляя (59), (60) в (58), прибавляя и вычитая в последнем div(XVT), получаем:
( 2 A ( 2 A ~
pu v рє 2 = -div рб u У w 2 -(6 ct')-XVT
V У V J J
+ PT+ 5Vsj - a'ik - divfcVT)
+
(61)
Сравнивая (61) с (53), получаем общее уравнение переноса тепла в вязкой жидкости:
pT(Is + '5VS) = + div(KVT) . (62)
Заметим, что при отсутствии вязкости и теплопроводности правая сторона (62) обращается в нуль и получается уравнение сохранения энтропии идеальной жидкости.
Общее уравнение переноса тепла (общее уравнение теплопроводности) (62) можно упростить в случае, когда скорость движения жидкости мала по сравнению со скоростью звука; возникающие в результате движения жидкости изменения давления при этом настолько малы, что вызываемыми ими изменениями плотности (и других термодинамических величин) можно пренебречь. Однако неравномерно нагретая жидкость не является все же при этом несжимаемой, так как ее плотность изменяется с изменением температуры. Поэтому при нахождении производных от термодинамических величин надо считать постоянным давление, а не плотность. Таким образом, имеем:
37
ds
dt
ds
dT
dT
dt ’
Vs =
ds
dT
VT
p
(63)
и ПОСКОЛЬКУ TI —
янном давлении, то
єсть теплоемкость при посто-
T f -- CP ї, ™ = CPVT ■ (64)
Уравнение (62) принимает вид
pCP+ BVT] = diV^VT^+&lk fX^ ■ (65)
Для того чтобы в уравнениях движения нагретой жидкости можно было считать плотность постоянной, необходимо (помимо малости отношения скорости жидкости к скорости звука), чтобы имеющиеся в жидкости абсолютные разности температур были достаточно малы. Тогда жидкость можно считать несжимаемой, т.е. применять уравнение непрерывности в виде divo = 0 ■
Считая абсолютные разности температуры малыми, можно пренебречь температурными изменениями величин 7, К , Cp , т.е. считать их постоянными. В результате можно получить уравнение переноса тепла в несжимаемой вязкой жидкости в сравнительно простом виде.
Анализируя закон сохранения энергии для несжимаемой жидкости, в которой имеют место процессы теплопроводности и внутреннего трения, можно относительно просто получить следующее уравнение переносатепла:
dT
dt
+ oVT = vAT h--
2Cp
'dol
ydXk
+
d»k )2
dXi J ,
(66)
n rs
где v = — — кинематическая вязкость; X ~—p---
p Р^-р
коэффициент температуропроводности; к — коэффициент теплопроводности; Cp — теплоемкость при постоянном давлении.
Если выделение тепла из-за вязкой диссипации незначительно, то (66) можно представить в виде
~RТ~ + ^VT = ZAT . (67)
dt
Особенно просто выглядит уравнение переноса тепла в неподвижной жидкости, где перенос энергии обязан целиком теплопроводности:
Ї (68)
Это выражение в математической физике называется уравнением теплопроводности. Оно может быть выведено и гораздо более простым образом, без помощи общего уравнения переноса тепла в движущейся жидкости.
6. Уравнения тепло- и массопереноса в жидкой смеси
Во всем предыдущем изложении предполагалось, что расплав (вязкая жидкость) однороден по своему составу. Если же расплав неоднороден, т.е. представляет собой некоторую смесь, состав которой
38
меняется вдоль ее объема, то уравнения гидродинамики существенно изменяются.
Ограничимся рассмотрением бинарной смеси (расплава, состоящего только из двух компонент). Бинарную смесь будем характеризовать концентрацией C - отношением массы одного из входящих в состав смеси веществ к полной массе смеси в данном элементе объема.
С течением времени распределение концентрации в расплаве меняется. Изменение концентрации происходит двумя путями. Во-первых, при макроскопическом движении расплава каждый данный элемент объема его передвигается как целое с неизменным составом. Этим путем осуществляется чисто механическое перемешивание расплава; хотя состав каждого перемещающегося элемента объема не меняется, но в каждой данной неподвижной точке пространства концентрация находящегося в этом месте расплава будет со временем изменяться. Если отвлечься от процессов теплопроводности и внутреннего трения, которые одновременно могут иметь место, то такое изменение концентрации является термодинамически обратимым процессом и не ведет к диссипации энергии.
Во-вторых, изменение состава может происходить путем молекулярного переноса веществ смеси из одного участка расплава в другой. Выравнивание концентрации путем такого непосредственного изменения состава каждого из участков расплава называют диффузией. Диффузия является процессом необратимым и представляет собой наряду с теплопроводностью и вязкостью один из источников диссипации энергии в жидкой смеси.
Будем обозначать посредством р полную плотность смеси. Уравнение непрерывности для полной массы смеси сохраняет прежний вид:
—— + divo = 0 dt
(69)
Оно означает, что полная масса жидкости в некотором объеме может измениться только путем втекания или вытекания жидкости из этого объема. Строго говоря, для жидкой смеси само понятие скорости должно быть определено заново. Уравнение непрерывности в виде (69) определяет скорость в соответствии с прежним определением как полный импульс единицы массы смеси.
Не меняется также и уравнение Навье-Стокса. Выведем остальные два гидродинамические уравнения для бинарной смеси.
В отсутствие диффузии состав каждого данного элемента смеси оставался бы неизменным при его передвижении. Это значит, что полная производ-
dC
ная -- была бы равна нулю, т.е. имело бы место
уравнение
dC
dt
dC
dt
+ oVC = 0
Это уравнение можно написать, используя (69), как
+ div(^^C) = 0 ,
РИ, 2001, № 4
т.е. в виде уравнения непрерывности для одного из веществ в смеси (рС есть масса одного из веществ смеси в единице обьема). Написанное в интегральном виде
^ j^CdV = -jpCudf,
V0 F0
оно означает, что изменение количества данного вещества в некотором обьеме V равно количеству этого вещества, переносимому движущейся жидкостью через поверхность обьема F0 .
При наличии диффузии наряду с потоком
jpCudf
данного вещества вместе со всей смесью имеется еще и другой поток, который приводит к переносу веществ в смеси даже при отсутствии движения жидкости в целом. Пусть j есть плотность этого диффузионного потока, т.е. количество рассматриваемого вещества, переносимого путем диффузии в единицу времени через единицу поверхности. Тогда для изменения количества этого вещества в некотором обьеме Vo имеем:
д_
dt
J pCdV
Vo
jpCddf - J Jdf
Fo
Fo
или в дифференциальном виде
^PC) _ -div^pCu) - divl dt
(70)
С помощью (69) это уравнение непрерывности для одного из веществ в смеси можно написать в виде
Р\
— + 5VC dt
= -divl
(71)
Для вывода еще одного уравнения необходимо повторить произведенный в п.5 вывод, учитывая, что термодинамические величины жидкости являются теперь и функциями концентрации. При
д_
вычислении производной
ри
+ ре
с помо-
2
щью уравнений движения в п.5 необходимо было
дє
преобразовать члены р— и (-5Vp). Эти преобра-
зования теперь изменяются в связи с тем, что термодинамические соотношения для энергии и тепловой функции содержат дополнительный член с дифференциалом концентрации:
Р
ds = Tds Л—— dp + LidC P2 ,
dw = Tds H— dp + pdC P
где p — химический потенциал смеси. Соответ-
дє
ственно этому в производную — войдет теперь
dC
дополнительный член . Написав второе из
термодинамических соотношений в виде
dp = pdw - pTds - ppdC,
видим, что в (- 5Vp) войдет дополнительное слагаемое ppoVC.
Поэтому к выражению (61) надо добавить
(dC _ б , 7
рр\—^ + uVC I = -pdivl .
В результате получим:
5 брб
dt { 2
+ рє I = -div
op
2
Л
■ + w
-(5 ст')-
(ds ^ Л dvt ^ - (72)
+ pTI----Ь uVs I - a'ik---Ь divq - pdivl
{dt ) dxk
Вместо (-XVT) теперь введен некий вектор плотности потока тепла q , который может зависеть не только от градиента температуры, но и от градиента концентрации. Сумму двух последних членов с правой стороны равенства представим в виде
divq - pdivl = div(q - pj^ +1Vр .
Выражение рб
^о2
■ + w
-(б ст')
+ q, стоящее под зна-
2
V /
ком дивергенции в (72), есть вектор плотности полного потока энергии в жидкости. Первый член характеризует обратимый поток энергии, связанный просто с перемещениемжидкости как целого, а сумма (-(бст')+ q) характеризует необратимый поток.
Уравнение закона сохранения энергии гласит:
д_
dt
( 2 > (2 >
po h рЄ 2 V У = -div рб О h W 2 V J -(б ст') + q
.(73)
Вычитая его почленно из (72), получаем некоторое уравнение
pT(Is + '5Vs) = a'k ~ diAq ~rf)~1, (74)
обобщающее выведенное ранее (62).
Таким образом, получена полная система гидродинамических уравнений для жидких смесей. Число их в этой системе на единицу больше, чем в случае чистой жидкости, соответственно тому, что имеется еще одна неизвестная функция-концентрация. Этими уравнениями являются: уравнение непрерывности (69), уравнение Навье-Стокса, уравнение непрерывности для одной из компонент смеси (70) и уравнение (74), определяющее изменение энтропии.
Далее следует перейти к нахождению векторов j и q. Диффузионный вектор плотности потока вещества j и тепловой вектор плотности потока q обусловлены наличием в жидкости градиентов концентрации и температуры. Если градиенты концентрации и температуры невелики, то можно считать, что j и q являются линейными функциями от Vр и VT . Соответственно этому имеем:
j =-aVp-fiVT , (75)
q = -5Vp - yVT + pj . (76)
РИ, 2001, № 4
39
Между коэффициентами р и S существует простое соотношение, являющееся следствием принципа симметрии кинетических коэффициентов [3]:
S = pT. (77)
Поэтому
I = -aVp~pVT , (78)
q = -pTVp-yVT + р/ . (79)
В выражении для q исключим градиент Vр, выразив его через / и ут :
q = | м+— \I-kvt
где
К = у-
р2т
(80)
(81)
Переходя к переменным p, T и C , найдем Уц :
Зц
ЗС
VC +
p,T
VT+М vp (82)
эт Jc,p Up )c,t ' <82)
Последнее слагаемое в (82) можно преобразовать, используя термодинамическое соотношение
dy = -sdT + Vdp + pdC , (83)
того, считать, что в жидкости нет никакого макроскопического движения помимо того, которое может быть вызвано самим наличием градиентов температуры и концентрации. Скорость такого движения будет пропорциональна этим градиентам, поэтому в уравнениях (71) и (74) члены, содержащие скорость, оказываются величинами второго порядка малости и могут быть опущены. Величиной второго порядка является также член /V р в (74). Таким образом, остается
р^С + div/ = 0 pT — + diviq - pf\ dt ’dt
Подставив сюда для / и q выражения (87) и (88) (без члена с Vp), получим после простых преобразований [3] следующие уравнения:
— = піде + AT
dt { T
(89)
І - СГ [§ ) ^дт
dt Cp \ dC ) pT
(90)
Эта система линейных уравнений определяет распределение температуры и концентрации в жидкости.
где у — термодинамический потенциал единицы массы; V — удельный объем. Имеем:
(dp\ д V Ґ дУ Л
{dp)CJ ~ dpdc ЧзсJp,T• (84)
Подставив Vp в (75) и введя обозначения
D =
a(dp\ P-tD ■ = a|
P\dC ) p,T , T •vt l. P (85)
(SV ) /( |
p p[acj pTI ^ ~dC , p,T , (86)
получим следующие выражения:
f k к
/ =-pD
VC + VT + '^ Vp
T p
\
q =
к і _Ф| _ т (9ц
T ’5е J pT T \дт Jp,С
+ ц
i -KVT
(87)
. (88)
Коэффициент D называется коэффициентом диффузии; он определяет диффузионный поток при наличии только градиента концентрации. Диффузионный поток, вызываемый градиентом температуры, определяется коэффициентом термодиффузии кtD ; безразмерную величину кт называют термодиффузионным отношением. В учете последнего слагаемого в (87) может возникнуть необходимость при наличии в жидкости существенного градиента давления, вызванного, например, внешним полем. Величину kpD можно назвать коэффициентом бародиффузии.
Далее ограничимся рассмотрением случая, когда нет существенного градиента давления, а концентрация и температура настолько мало меняются в жидкости, что коэффициенты в выражениях (87), (88), являющиеся в общем случае функциями от С и т , можно считать постоянными. Будем, кроме
40
В особенности важен случай, когда концентрация смеси мала. При стремлении концентрации к нулю коэффициент диффузии стремится к некоторой конечной постоянной, а коэффициент термодиффузии — к нулю. Поэтому при малых концентрациях кт мало, и в уравнении (89) можно пренебречь членом ктVT . Оно переходит тогда в уравнение диффузии:
f=DAC- (91)
7. Уравнение Навье-Стокса и уравнение непрерывности в переменных вихрь, функция тока
Для однородной по своему составу несжимаемой вязкой жидкости исходная система уравнений состоит из уравнения Навье-Стокса и уравнения непрерывности:
4^ + (бУ)б = --^gradp + vAu , (92)
div о = 0 . (93)
В этой системе искомыми функциями являются вектор скорости б и давление p, которые зависят от пространственных координат и времени t, а параметрами будут плотность р и кинематическая вязкость v.
Для класса двумерных движений в декартовой системе координат x, y система уравнений (92), (93) принимает вид
dux dux dux 1 dp
x- + »х~:г + »y — =-~ + v
f P2
dt
dx
dy p dx
,2 A
d ox ^ d ox
doy doy doy
—-+°x^:r+uy —
dt
dx
1 dp ■-----— + v
V ( V2
dx2 су2
dy p dy
do
J
У+д 4 ^
22
ox dy
V J J
,(94)
,(95)
РИ, 2001, № 4
дох
дх
= 0
(96)
При выполнении вычислений обычно используется безразмерная форма записи исходных уравнений. Для этого вводят безразмерную величину в виде
(р =
П
где (pi — некоторый масштаб.
(97)
Безразмерная запись уравнений может быть получена подстановкой в уравнение вместо размерной величины (р ее выражения в виде ycpj. Выбор масштабов зависит от конкретной постановки задачи. Пусть, например, условиями задачи заданы характерная скорость Oj и размер области L . Течение с заданной характерной скоростью или перепадом давления называется вынужденной конвекцией. Выбирая в качестве масштабов скорости и длины соответственно oj, L , а для параметров р , v — их значения, заданные условиями задачи, получаем следующую безразмерную запись исходной системы уравнений:
+ (иV)б = -gradp + R До, (98)
divd = 0 . (99)
В этой системе имеется безразмерный параметр, построенный по величинам, выбранным в качестве масштабов: Re = UjL lv — число Рейнольдса, представляющее отношение сил инерции к силам вязкости и определяющее интенсивность вынужденной конвекции.
Заметим, что для времени ~ и давления p в системе (98), (99) использованы масштабы в виде
L
^ и ро2 соответственно.
Использование безразмерной системы преследует две цели: приведение значений вычисленных величин к соответствующей шкале, а также расчет и обработка результатов в форме, содержащей минимальное число параметров.
Уравнения (94)-(96) могут быть записаны в иной форме, не содержащей давления и в ряде случаев более удобной для численной реализации. Из уравнений (94), (95) можно исключить давление p, продифференцировав первое из них по у , а второе по x . Определяя вихрь как
дОх доу
ду дх ,
(100)
получаем уравнение переноса вихря:
да
dt
~ + Ох
да
дх
- + оу
да
ду
<д 2а
— V
д2^
бх2 ду1
(101)
Следует заметить, что в трехмерном случае вихрь обычно определяется как [v б], что при переходе к двухмерному случаю дает выражение, отличающееся от (100) знаком.
Определяя функцию тока у соотношениями
ду
ду
= -о
у J
(102)
9у ’ дх
получаемуравнение связи вихря а и функции тока у:
д 2у д 2у
Р? + у (103>
Отметим, что структура соотношений (102) следует из уравнения непрерывности для несжимаемой жидкости:
, ^ бо х diw =
дх
5о у
+ —21 = о бу '
(104)
Зная функцию тока, можно определить форму его линий. Действительно, дифференциальное уравнение линий тока (при двумерном течении) есть
йх dy
vx =Vy (105)
или оуйх -охйу = 0 ; оно выражает тот факт, что направление касательной к линии тока в каждой точке совпадает с направлением скорости. Подставляя сюда (102), получаем:
ду ду
—— йх Л—— йу = dy = 0
дх ду '■
(106)
откуда у = const. (107)
Таким образом, линии тока представляют собой семейство кривых, получающихся приравниванием функции тока у{х, у) произвольной постоянной.
Итак, в переменных Д, у система уравнений (94)-(96) принимает вид
да
dt
■ + Ох
да
дх
- + оу
да
ду
fd^
дх 2
sV
бу 2
5 2у д2у -+—Ь- = а
дх2 ду2
8. Единая математическая модель неоднородной жидкости в приближении Буссинеска
Плотность и физические свойства неоднородной жидкости изменяются в зависимости от пространственной переменной. Причиной неоднородности жидкости может быть изменение ее состава или температуры, что приводит к ряду новых эффектов, которые отсутствуют в однородной изотермической жидкости (конвекция, тепло- и массоперенос).
Рассмотрим бинарную смесь, уравнение состояния которой задано в виде
р = Лс, т), (110)
где т — температура; C — концентрация примеси; предполагается, что плотность р не зависит от давления p.
Далее предположим, что
Р = Р0 + p’, Р = Р0 +р’, (111)
2
(108)
(109)
где Р0, р0 удовлетворяют уравнению гидростатики:
gradP0 = Р0g , (112)
Р'<< Р0 , р’«Р0. (113)
РИ, 2001, № 4
41
Тогда, с точностью до малых первого порядка имеем:
Yl+v_J_ vpoP’
Р А) А) А)
С учетом (112), перепишем (114) в виде
Vp Vp _ 1 ,
---=----+ g---2 Vp0 p
P Po PQ
(114)
(115)
Преобразуем третье слагаемое в правой части (115). Напишем переменные т и C в виде
т = т0 + T', (116)
C = C0 + C', (117)
где Tо , Cо — некоторые постоянные средние значения, от которых отсчитываются неравномерности температуры T' и концентрации C’, причем
т' << т0 , C' << C0 .
Ввиду малости изменения температуры т' и концентрации C' мало также изменение плотности р', поэтому можно написать:
Р =
дрр
дт
т' +
p,C
дР0
dC
C'
р,т
и соответственно
(118)
—у VP0P'
Р 0
_g_
Р0
= Рт^' + PcgC' >
_Ф^| т, g_ Гбр) C ’ =
дт У 0; p,C р0 UC у 0; р,т
(119)
здесь Рт =—— (120)
Р0 У^т)0;p C
— коэффициент теплового изменения плотности,
PC dC] 0;p, C “21)
— коэффициент концентрационного изменения плотности.
Теперь перейдем к преобразованию уравнения НС при наличии поля тяжести:
до /-„w Vp - -
— + \uVp = —— + vAu + g . (122)
Поскольку
g + Ртііт' + PcgC', (123)
Р A) v '
то уравнение Н-С (122) принимает (при опущенном индексе у р0) вид:
Vp' + vAo + ^Ртт' + РсС 0. (124)
Уравнение движения (124) совместно с уравнением непрерывности, уравнением переноса тепла и массы неоднородной жидкости составляют систему уравнений, представляющую математическую модель неоднородной жидкости в приближении Бус-синеска:
divo = 0 , (125)
pCp
дт ’ ~dt
(oV )т'
= кдт'
(126)
дС'
—+ (6V)C' = DAC’. (127)
dt
В системе уравнений (124) - (127) искомыми являются вектор скорости (3 , давление p , температура т', концентрация примеси С' (точнее отклонения последних от их статических значений), которые зависят от пространственных координат и времени t. Параметрами являются: плотность р, кинематическая вязкость v, коэффициенты теплопроводности К и диффузии D, удельная теплоемкость при постоянном давлении cp , ускорение свободного падения g и коэффициенты теплового и концентрационного изменения плотности Рт и Рс .
Для решения системы (124)-(127) необходимо задать начальные и граничные условия.
Следует отметить так называемые граничные условия прилипания, в соответствии с которыми все компоненты скорости на твердой неподвижной стенке равны нулю: ог = 0 .
Граничные условия для полей температуры и концентрации могут быть трех основных типов:
1) задана температура Тг (концентрация примеси Сг) на границе плоской поверхности жидкости;
_ ^дт
2) задан поток тепла qr --N—— (диффузионный
, „дС дп
поток I г --D——), n — нормаль к поверхности
on
жидкости;
3) задан закон теплообмена в виде qr =ат (тг - тс), где ат — коэффициент теплопередачи; Тг — температура стенки; тс — температура среды (соответственно закон массообмена в виде Iг = ас(Сг - Сс), где ас — коэффициент массообмена).
Вводя масштабы для искомых величин и независимых переменных, можно привести систему к безразмерному виду.
Относительно просто представима система (124) -(127) в прямоугольной и цилиндрической системах координат в естественных переменных и переменных СО , V .
9. Заключение
Детально исследована единая математическая модель процесса выращивания монокристаллов методом Чохральского. Конкретные формулировки задач выращивания монокристаллов основаны на этой модели и соответствующих начальных и граничных условиях.
Литература: 1. Лодиз Р., Паркер Р. Рост монокристаллов. М.: Мир, 1974. 540 с. 2. Шашков Ю.М. Выращивание монокристаллов методом вытягивания. М.: Металлургия, 1982. 312 с. 3. Ландау Л.Д., Лифшиц Е.М. Гидродинамика. М.: Наука, 1986. 736 с. 4. Математическое моделирование. Получение монокристаллов и полупроводниковых структур / А.А. Самарский, Ю.П. Попов и др. М.: Наука, 1986. 200 с. 5. Математическое моделирование конвективного тепломассообмена на
РИ, 2001, № 4
42
основе уравнений Навье-Стокса / В.И. Полежаев, А.В. Бунэ, Н.А. Верезуб и др. М.: Наука, 1987. 272 с.
Поступила в редколлегию 11.07.2001 Рецензент: д-р физ.-мат. наук Гордиенко Ю.Е.
Тевяшев Андрей Дмитриевич, д-р техн. наук, профессор, зав. каф. ПМ ХНУРЭ. Научные интересы: системный анализ и теория оптимального стохастического управления. Адрес: Украина, 61166, Харьков, пр. Ленина 14.
Суздаль Виктор Семенович, канд. техн. наук, зав. отделом НТК “Институт монокристаллов”. Научные ин-
тересы: системы управления технологическими процессами получения монокристаллов. Адрес: Украина, 61001, Харьков, пр. Ленина, 60.
Бородавко Юрий Михайлович, канд. техн. наук, профессор кафедры ПМ ХНУРЭ. Научные интересы: математическое моделирование физических процессов. Адрес: 61166, Харьков, пр. Ленина, 14.
Пелипец Андрей Александрович, аспирант кафедры ПМ ХНУРЭ. Научные интересы: математическое моделирование процессов выращивания монокристаллов из расплава. Адрес: Украина, 61166, Харьков, пр. Ленина 14.
УДК 548.771.319
МЕТОД КОНТРОЛЯ УРОВНЯ РАСПЛАВА ПРИ ВЫРАЩИВАНИИ СТРУКТУРНО-СОВЕРШЕННЫХ МОНОКРИСТАЛЛОВ КРЕМНИЯ ПО МЕТОДУ ЧОХРАЛ ЬСКОГО
ОКСАНИЧ А.П, ПРИТЧИН С.Э., ПЕРВАКОВА.С.
Из всех разработанных на сегодня методов контроля уровня расплава при выращивании структурно-совершенных кристаллов кремния наиболее перспективен метод анализа телеизображения поверхности расплава вблизи зоны роста монокристалла. Он обеспечивает необходимую точность измерений при простоте конструкции измеряющего устройства и широких возможностях интеграции в комплекс устройств контроля за условиями выращивания монокристаллов.
1. Введение
В процессе выращивания монокристаллического слитка из расплава по методу Чохральского изменяется уровень расплава вблизи области, где происходит рост монокристалла. Основными причинами этого являются:
1. Снижение уровня расплава вследствие переноса массы на выращенный монокристалл.
2. Колебания зеркала расплава, вызванные столкновением с ним охлаждающих потоков аргона. Используемые для выращивания монокристалла по методу Чохральского установки обычно имеют ростовую камеру, в которой во время процесса сверху вниз движется поток аргона, охлаждающий выращиваемый кристалл. Часть потока аргона, однако, сталкивается с поверхностью расплава, вызывая тем самым колебания его зеркала относительно среднего уровня.
Необходимость контроля уровня расплава обусловлена следующими причинами:
1. Колебания поверхности расплава усложняют продолжение выращивания кристалла (особенно во время роста шеечной части, в которой диаметр кристалла мал) так как нижняя часть выращенного кристалла оказывается под угрозой отрыва от поверхности расплава, что сделает невозможным дальнейшее выращивание кристалла. На участке выращивания шейки, после погружения в расплав затравки,
РИ, 2001, № 4
до формирования конического участка требуется предельная осторожность, так как кристалл выращивается со сравнительно высокой скоростью порядка 2 мм/мин или более, когда его диаметр ограничен 2-5 мм. На практике используется ручной контроль, что приводит к тому, что около 10% таких результатов оказываются неудачными даже в случае, если их осуществляют опытные специалисты. Когда диаметр выращиваемого кристалла слишком умень -шается, его нижняя часть отрывается от поверхности расплава, делая невозможным последующее его выращивание.
2. Для получения бездефектных монокристаллических слитков кремния необходимо поддерживать постоянство регулирующих параметров, таких как скорость роста и температура расплава, в течение длительного периода времени. Для удовлетворения этих требований необходима неизменность уровня расплава относительно нагревательного узла ростовой камеры.
На зависимость отводимой из системы кристаллрасплав теплоты кристаллизации указано в [1]. Тепловой баланс можно представить в следующем виде:
dm
dT
L------ъ k-----A- — k
dt
dx-i
dT
dx2
A-)
(1)
где L — скрытая теплота кристаллизации;
dm
dt
доля закристаллизовавшегося в единицу времени расплава; T—температура расплава; Гд— теплопро-
dT
водность расплава;
dx
- температурный градиент
в расплаве в некоторой точке х1 вблизи поверхности раздела; Ai — площадь изотермической поверхности, проходящей через х1; ks — теплопроводность
dT
кристалла; dx — температурный градиент в слитке
в точке х2 вблизи поверхности раздела; А2 — площадь изотермической поверхности, проходящей через х2, примерно равная сечению кристалла; если изотермическая поверхность плоская и перпендикулярна к направлению роста, она точно соответствует сечению.
При этом:
dm. = vAc8 = Vp + VdVc , (2)
здесь v — скорость роста, равная сумме скорости вытягивания Vp и скорости опускания поверхности
43
