Investigation of Marangoni convection during contactless crystal growth in microgravity conditions

Толық мәтін

ВВЕДЕНИЕ

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

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

Ранее нами были опубликованы результаты исследований реальной структуры кристалла GaSb:Te, выращенного на борту автоматического космического аппарата методом охлаждения в градиенте температуры [5–7]. Внешние размеры ампулы составляли: диаметр 10 мм, длина 65 мм, размеры образца — диаметр 6 мм, длина 30 мм. Исходная концентрация Te в кристалле составляла 5.6 × 10¹⁷ см^–3. Рост проводился в направлении [111]B. Примерно половина кристалла выросла без контакта со стенками ампулы и обладала весьма высоким структурным совершенством [8]. В [9] методом количественной плосковолновой рентгеновской топографии [10, 11] была построена карта вариаций концентрации Te на поверхности продольного среза кристалла, которая продемонстрировала хорошее согласие с данными фотолюминесценции [12]. Это дает возможность количественного анализа влияния конвекции на захват примеси растущим кристаллом.

В [13] было показано, что на основе одномерных аналитических моделей, описывающих сегрегацию примесей в кристалле в процессе роста возможно определить параметры ростового процесса, например, скорость роста и скорость конвекции, по данным о распределении примеси в кристалле. Известны стационарные решения в моделях Бартона–Прима–Слихтера (БПС) [14], Острогорского–Мюллера (OM) [15] и Garandet [16]. Полные решения даны в работах [13] (модели БПС и OM) и [17] (модель Garandet).

Точность моделей БПС, OM и Garandet исследовалась в [18, 19] в сопоставлении с результатами двумерного численного моделирования. Было показано, что погрешность всех моделей находится в интервале ±10%, однако наиболее надежной и точной оказалась модель БПС, погрешность которой не превысила 5%. Таким образом, эти простые модели оказываются вполне пригодными для количественных оценок параметров роста кристалла, в том числе скорости конвекции, по данным о содержании примеси в кристалле.

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

АНАЛИЗ РАСПРЕДЕЛЕНИЯ ПРИМЕСИ В КРИСТАЛЛЕ

Некоторые свойства кристалла GaSb и его расплава приведены в табл. 1. Значения удельной свободной энергии поверхности раздела расплав–кристалл σ в литературе найти не удалось, и мы попытались оценить его, опираясь на известный факт, что для алмазоподобных полупроводников существует линейная зависимость величин σ и σ_s от среднего атомного номера соединения [34]. Используя значения σ = 340 эрг/см для Si [35] и σ = 200 эрг/см для Ge [36], методом экстраполяции мы определили для GaSb σ ≈130 эрг/см. При использовании значения k_n0 = 0.37 для примеси Te в кристалле GaSb получаем k_f0 = 3.5, то есть k_n0/k_f0 ≈ 9.5.

 

Таблица 1. Свойства кристалла и расплава GaSb

Параметр	Обозначение	Значение	Ссылка
Температура плавления,°C	Tm	712	[20]
Угол роста	e	10° 15° 28°	[21] [22] [23]
Угол контакта расплава с поверхностью SiO2 (кварц)	am	121°	[24]
Угол контакта расплава с поверхностью C (графит)	aС	128°	[24]
Поверхностное натяжение на поверхности раздела расплав–газ при 725°C (эрг/см)	sm	454 453 450	[25] [24] [23]
Удельная энергия поверхности раздела кристалл–газ для грани (111) (эрг/см)	sS	910	[26]
Удельная энергия поверхности раздела расплав–кристалл (эрг/см)	s	~130	Оценено по среднему атомному номеру
Температурная зависимость sm (эрг/см ⋅К)		0.11 0.098±0.027	[25]
Плотность расплава (кг/м3) при 712°C при 1000°C	r	6.05 × 103 5.87 × 103	[12]
Коэффициент диффузии Te в кристалле при Tm (см2/с)	DS	1.34 × 10–9	Рассчитано по [27]1
Коэффициент диффузии Te в расплаве при Tm (см2/с)	Dl	3 × 10–5 5 × 10–5	[28] [29], [20]
Теплоемкость кристалла при 298 K (эрг/г⋅K)	Cps	2.54 × 106	[20]
Теплоемкость расплава при Tm (эрг/г⋅K)	Cpl	3.998 × 106	[30]
Коэффициент теплопроводности расплава при Tm (эрг/см⋅с⋅K)	cl	21.7 × 105	[31]
Коэффициент теплопроводности кристалла при Tm (эрг/см⋅с⋅K)	cS	7.81 × 105	[31]
Коэффициент температуропроводности расплава при Tm (см2/с)	ll	9 × 10–2	—
Коэффициент температуропроводности кристалла при Tm (см2/с)	lS	5.5 × 10–2	—
Кинематическая вязкость расплава при Tm (см2/с)	n	3.75 × 10–3	[30]
Динамическая вязкость расплава при Tm (г/см⋅с)	m	2.261 × 10–2	[30]
Объемный коэффициент температурного расширения расплава (K–1)	bT	9.6 × 10–5	[32]
Кинетический коэффициент ступени на грани (111) (см/с⋅К)	b	18	[33]

¹ В [28] при вычислении диффузионных параметров Te в кристалле GaSb (предэкспоненциального множителя D₀ и энергии активации Q) была допущена ошибка, которая затем была перенесена в книгу “Атомная диффузия в полупроводниках”, п/ред Д.Шоу. М.: Мир. 1975. 684 С. Правильные значения диффузионных параметров следующие: D₀ = 1.53 × 10^–2 см²/с, Q = 2.21 × 10^–12 эрг.

 

На рис. 1 представлена карта вариаций концентрации Te, полученная методом плосковолновой рентгеновской топографии. Область A образована ростом округлого фронта, а область B — плоской гранью. Поскольку в области роста грани коэффициент захвата примеси отличается от равновесного, количественный анализ проводился по области A с использованием данных [9].

 

Рис. 1. Области кристаллизации округлого фронта (А) и грани (В).

 

В [19] показано, что одномерные аналитические модели адекватно описывают распределение примеси в кристалле, если значение ее концентрации в точке принимается равным усредненной по сечению концентрации вдоль соответствующего перпендикулярного среза кристалла. Для округлого фронта кристаллизации усреднение концентрации Te проводилось по области A (рис. 1). Результаты представлены на рис. 2.

 

Рис. 2. Усредненное параллельно фронту кристаллизации распределение Te в области роста округлого фронта и расчет по модели БПС. В верхней части рисунка приведен график изменения зазора δl между стенкой ампулы и левой стороной кристалла со сдвигом на 0.5 мм (пояснения см. в тексте).

 

При анализе исходная концентрация примеси в расплаве c₀ принималась равной средней концентрации Te в затравке (нерасплавленной части исходного монокристаллического слитка). Пусть x — это координата в неподвижной системе координат, обозначающая расстояние от точки начала кристаллизации, и вычисляемая как

$x=Rt+\xi ,$ (1)

где R — нормальная скорость роста кристалла, t — время, ξ — координата в подвижной системе координат, связанной с фронтом кристаллизации, обозначающая расстояние от фронта кристаллизации.

На рис. 2 видно, что начальный участок роста кристалла (до x ≈ 0.48 см) хорошо описывается в рамках модели БПС:

$k\left(x\right)=\frac{c_S\left(x\right)}{c_0}=k^{*}-(k^{*}-k_0)e^{-\alpha x}$, (2)

где

$\alpha =\frac{R}{D}\frac{k_0}{k^{*}}\frac{1}{1-(1+\Delta )e^{-\Delta }}$, (3)

$k^{*}=\frac{k_0}{k_0+(1-k_0)e^{-\Delta }},\Delta =\frac{R\delta }{D}$, (4)

$\delta =\delta {}_{\nu }(D/\nu )^{1/2}=4.6{(\nu L/V_{\infty })}^{1/2}{(D/\nu )}^n$ [37], (5)

со следующими параметрами:

– начальная концентрация примеси в жидкой фазе c₀ = 5.6 × 10¹⁷ см^–3,

– концентрация примеси в кристалле в стационарном режиме $c_S^{*}$ = 5.15 × 10¹⁷ см^–3,

– эффективный коэффициент распределения примеси в стационарном режиме k* = 0.92,

– R = 4.2 × 10^–4 см/с,

– скорость конвективного потока за пределами пограничного слоя V_∞ = 8.4 × 10^–3 см/с.

Поскольку одномерные аналитические модели построены в предположении постоянства скоростей роста кристалла и конвективных потоков, то, строго говоря, они не могут быть применены для анализа участка кристалла при x > 0.48 см, где наблюдаются значительные вариации эффективного коэффициента распределения примеси в стационарном режиме k*. Однако, считая, что к данному моменту рост кристалла вышел на стационарный режим, а последующие изменения k* происходили достаточно медленно, можно в первом приближении использовать стационарное решение в модели БПС для оценки изменения параметров процесса. Как уже указывалось, такое решение зависит от отношения скоростей конвекции и роста кристалла и не позволяет определить величины этих параметров одновременно. Поэтому для их оценки мы будем привлекать дополнительные соображения.

Уровень остаточной микрогравитации для различных космических аппаратов составляет величину ~10^–5–10^–4 g/g₀ (g₀ = 980 см/с), что в нашем случае должно вызывать конвекцию со скоростями порядка ~10^–6–10^–4 см/с [18]. Полученное нами значение скорости V_∞ = 8.4 × 10^–3 см/с минимум на полтора порядка больше, что заставляет предположить превосходящий вклад термокапиллярной конвекции (конвекции Марангони), обусловленной градиентом поверхностного натяжения расплава в условиях неоднородного теплового поля.

В условиях бесконтактного роста кристалла, когда его диаметр меньше внутреннего диаметра ростового контейнера, вблизи фронта кристаллизации образуется участок свободной поверхности расплава, называемый мениском (основная часть расплава находится в контакте со стенками контейнера). В процессе роста кристалла мениск движется вместе с фронтом кристаллизации. Такая модель была впервые предложена в [38, 39] и впоследствии обсуждалась в [40–42] (главным образом, в плане происхождения мениска).

Высота мениска (рис. 3) зависит от его формы, расстояния между кристаллом и стенкой контейнера δl, угла смачивания расплавом стенки контейнера a_m и угла роста e = j — q (q — угол между поверхностью фронта кристаллизации и боковой поверхностью кристалла, j — угол между поверхностью мениска и фронтом кристаллизации). Точка О – это точка пересечения трехфазной линии с плоскостью рисунка. Значения углов a_m и e вместе с другими параметрами для GaSb приведены в табл. 1.

 

Рис. 3. Геометрия межфазных границ вблизи фронта кристаллизации при отсутствии контакта кристалла со стенками контейнера.

 

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

$h_m=\delta l\frac{\sin {\alpha }_m+\cos \phi }{\cos {\alpha }_m+\sin \phi }$ . (6)

Для случая роста кристалла постоянного диаметра (q = 90°, j = 90° + e) имеем

^{$h_m=\delta l\cdot tg\frac{{\alpha }_m-\epsilon }{2}$.} (7)

Между левой стороной образца и стенкой ампулы имелся существенный зазор, размер которого менялся со временем. Сопоставление графиков (рис. 2) показывает, что между изменением величины зазора dl и концентрацией Te в образце имеется достаточно отчетливая корреляция.

Распространение движения жидкости с поверхности расплава внутрь его объема связано с действием сил вязкого трения и развитие этого процесса требует времени ~ L²/v (L — характерный размер). Соответственно, между изменением величины зазора и соответствующим изменением содержания Te в образце в нашем случае (при L = 0.6 см) должно пройти время ~100 с, что при скорости роста R ≈ 4–5 × 10^–4 см/с составит длину прироста кристалла ~ 0.4–0.5 мм. График изменения dl (рис. 2) дан со сдвигом на 0.5 мм, при котором на наш взгляд, обсуждаемые ниже особенности процесса проявляются наиболее отчетливо. Заметим, что обсуждаемая ниже картина сохраняет все принципиальные черты в довольно широком интервале смещений графика dl относительно начала координат — от 0.2 до 1.0 мм. При этом меняются (в пределах 30 %) лишь получаемые величины скоростей роста отдельных участков кристалла и соотношение между ними.

Итак, для точек при x > 0.48 см применим формулу (4), из которой находим:

$\Delta \left(x\right)=-\ln \left(\frac{k(1-k^{*}(x\left)\right)}{k^{*}\left(x\right)(1-k)}\right),\text{\hspace{1em}где\hspace{0.33em}}{k}^{*}\left(x\right)=\frac{c_{S\left(x\right)}}{c_0}$, (8)

где c₀ = 5.6 × 10¹⁷ см^–3. График зависимости $\Delta \left(x\right)$ от dl представлен на рис. 4. Наиболее явная корреляция между вариациями концентрации примеси в кристалле и величиной dl наблюдается на участках 0.48 < x < 0.7 см и 0.75 < x < 0.85 см (рис. 2). Соответствующие этим участкам значения $\Delta \left(x\right)$ показаны на рис. 4 черным цветом. Эти точки хорошо аппроксимируются зависимостью:

${\Delta }_0\left(x\right)=\delta \left(x\right)\frac{R_0}{D}=\frac{\mathrm{0,047}}{\delta l\left(x\right)}$, (9)

где через ${\Delta }_0$ обозначена величина параметра $\Delta$в модели БПС при постоянной скорости R = R₀. Будем полагать, что на этих участках вариации концентрации примеси обусловлены вариациями скорости конвективного потока, а скорость роста кристалла не меняется и равна R₀ = 4.2 × 10^–4 см/с.

 

Рис. 4. Зависимость $\Delta \left(x\right)$ от δl, рассчитанная по измеренным значениям C_Te , и ее аппроксимация по формуле (9).

 

Для тех участков кристалла, где вариации концентрации Te не коррелируют с вариациями dl (рис. 2), значения $\Delta \left(x\right)$ на рис. 4 представлены красным цветом. Поскольку $\Delta \left(x\right)$ зависит только от R и V_∞, а величину V_∞ мы связываем с величиной dl, можно полагать, что участки , отклоняющиеся от зависимости (9), обусловлены вариациями скорости роста кристалла R.

С учетом (5) и (9) получаем

$V_{\infty }\left(x\right)={\left(98R_0\delta l\right(x\left)\right)}^2\frac{L}{D}$. (10)

Поскольку из (9) вытекает, что $\frac{\delta \left(x\right)}{D}=\frac{\mathrm{0,047}}{R_0\delta l\left(x\right)}$, а $R\left(x\right)=\Delta \left(x\right)\frac{D}{\delta \left(x\right)}$, то

$R\left(x\right)=21\Delta \left(x\right)\delta l\left(x\right)R_0$, (11)

где R₀ = 4.2 × 10^–4 см/с, L = 0.6 см.

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ЗАХВАТА ПРИМЕСИ КРИСТАЛЛОМ ПРИ НАЛИЧИИ КОНВЕКЦИИ В РАСПЛАВЕ

Для проверки найденных зависимостей было проведено 2D-моделирование движения расплава GaSb в условиях наличия его свободной поверхности вблизи фронта кристаллизации.

Рассматривалась модель конвективного тепломассопереноса в расплаве, предполагающая задание плоского движущегося фронта кристаллизации с учетом потока примеси из расплава в кристалл. Моделировалась кристаллизация с постоянной скоростью R = 3 × 10^–4 см/с плоского слоя расплава GaSb:Te высотой L и длиной H, находящегося в гравитационном поле g при T₀ = 985 K (температура плавления) и различных величинах T_w = 996–1057 K (рис. 5), соответствующих вариациям продольного температурного градиента.

 

Рис. 5. Схема, иллюстрирующая модель для численного решения двумерной задачи конвективной диффузии.

 

Векторная форма уравнений Навье–Стокса–Буссинеска переноса тепла и примеси в расплаве записывается в системе координат, связанной c движущимся фронтом кристаллизации, и имеет следующий вид [43]:

∂V/∂t + ((V–R)∇)V = –1/r∇P + nDV + gb_TT, (12)

div V = 0, (13)

rc_pl(∂T/∂t + ((V–R)∇)T) = lDT, (14)

∂c_L /∂t + ((V–R)∇)c_L = DDc_L. (15)

В этой системе уравнений искомыми являются вектор скорости потока расплава V, вектор скорости кристаллизации R = (R,0), давление P, температура T, концентрация примеси c_L, которые зависят от пространственных координат и времени t.

Для уравнений (12) — (15) задавались следующие начальное и граничные условия:

на границах расплава V = 0 и T задается в соответствии с рис. 5, (16)

начальное распределение примеси:

c_L = c₀ при t = 0, (17)

на фронте кристаллизации (x = 0):

$D\nabla c_L= \left(1-k\right)c_L$, (18)

на всех остальных границах:

$\nabla c_L= 0$. (19)

Итерационно решалась двумерная задача (12) — (19) в прямоугольнике x = (0, H), y = (0, L), в расчетах использовалась явная разностная схема 2-го порядка точности.

Расчеты для ячейки размерами L = 0.6 см, H = 3.0 см проводили при нулевом уровне гравитации, скорости роста кристалла R = 3 × 10^–4 см/с и осевом градиенте температуры G = 3.7 K/см. Величину δl принимали равной 0.01, 0.02, 0.03 и 0.05 см, мениск моделировали двумя симметричными участками боковой поверхности расплава высотой h_m (рис. 6), которая определялась по формуле (7) при значениях углов α_m = 121°, e = 28° (табл. 1).

 

Рис. 6. Схема, иллюстрирующая модель для численного расчета скорости конвекции Марангони.

 

Результаты расчетов представлены на рис. 7. Полученные данные неплохо аппроксимируются квадратичной зависимостью $V_{\infty }\left(x\right)=20\delta l{\left(x\right)}^2$ (рис.7а). Подстановка параметров, при которых проводилось моделирование, в формулу (10) дает близкую зависимость $V_{\infty }\left(x\right)=17\delta l{\left(x\right)}^2$.

 

Рис. 7. Результаты 2D моделирования конвекции Марангони: а – зависимость максимальной скорости конвективного потока V_∞ от зазора δl между боковой поверхностью кристалла и стенкой контейнера; б — распределение скорости конвективного течения; в – распределение температуры при δl = 0.05 см.

 

Таким образом, проведенные выше оценки представляются вполне разумными, что позволяет использовать формулы (10) и (11) для расчётов $V_{\infty }\left(x\right)$ и $R\left(x\right)$.

ЗАКЛЮЧЕНИЕ

Анализ начального переходного режима роста кристалла в рамках модели БПС позволил определить параметры роста кристалла на начальном (до выхода на стационарный режим) участке роста кристалла: c₀ = 5.6 × 10¹⁷ см^–3, $c_S^{*}$ = 5.15 × 10¹⁷ см^–3, k* = 0.92, R = 4.2 × 10^–4 см/с, V_∞ = 8.4 × 10^–3 см/с. Последующие вариации эффективного коэффициента распределения могли быть вызваны колебаниями как скорости роста кристалла, так и скорости конвективного потока. Поскольку стационарное решение не позволяет определить оба этих параметра одновременно, дальнейший анализ проводился в предположении об определяющем вкладе конвекции Марангони в скорость потока расплава. В рамках этой гипотезы скорость конвективного потока определяли по величине зазора между боковой поверхностью кристалла и стенкой ростового контейнера на основе эмпирической зависимости. Это позволило определить как величину V_∞, которая достигала значения 6.0 × 10^–2 см/с, так и скорость роста кристалла R_n в области округлого фронта. Было установлено, что, начиная с некоторого момента, R_n стала возрастать (до значения ≈8 × 10^–4 см/с), но не монотонно — области относительно плавного возрастания R_n прерывались короткими участками сравнительно резкого ее снижения.

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

ФИНАНСИРОВАНИЕ РАБОТЫ

Работа выполнена в рамках выполнения государственного задания НИЦ “Курчатовский институт” и в части численного моделирования по теме государственного задания ИПМ РАН № 124013000674-0.

КОНФЛИКТ ИНТЕРЕСОВ

Авторы данной работы заявляют, что у них нет конфликта интересов.

Автор А.Э. Волошин является членом редколлегии журнала “Поверхность. Рентгеновские синхротронные и нейтронные исследования”.

Авторлар туралы

А. Voloshin

National Research Centre “Kurchatov Institute”; Mendeleev Russian University of Chemical Technology; National University of Science and Technology “MISIS”

Хат алмасуға жауапты Автор.
Email: voloshin@crys.ras.ru

Shubnikov Institute of Crystallography, Kurchatov Complex of Crystallography and Photonics

Ресей, Moscow; Moscow; Moscow

Е. Rudneva

National Research Centre “Kurchatov Institute”

Email: voloshin@crys.ras.ru

Shubnikov Institute of Crystallography, Kurchatov Complex of Crystallography and Photonics

Ресей, Moscow

V. Manomenova

National Research Centre “Kurchatov Institute”

Email: voloshin@crys.ras.ru

Shubnikov Institute of Crystallography, Kurchatov Complex of Crystallography and Photonics

Ресей, Moscow

А. Prostomolotov

Ishlinsky Institute for Problems in Mechanics of Russian Academy of Sciences

Email: voloshin@crys.ras.ru
Ресей, Moscow

N. Verezub

Ishlinsky Institute for Problems in Mechanics of Russian Academy of Sciences

Email: voloshin@crys.ras.ru
Ресей, Moscow

Әдебиет тізімі

Қосымша файлдар