Содержание к диссертации
Введение
Алгоритм численного решения трехмерных уравнений Навье- Стокса для несжимаемой жидкости в естественных переменных 26
1.1 Система уравнений 26
1.2 Расщепление по физическим процессам 28
1.3 Численный метод решения уравнений Навье-Стокса для несжимаемой жидкости 29
1.4 Анализ различных способов решения уравнения для давления 34
1.5 Тестовые расчеты 38
1.6 Основные результаты 47
Конвективная устойчивость 48
2.1 Основные предположения 48
2.2 Математическая модель 50
2.3 Метод решения задачи концентрационной конвекции . 52
2.4 Модельная задача 54
2.5 Некоторые дополнительные теоретические и экспериментальные данные 60
2.6 Численное исследование возникновения подкритической конвекции 61
2.7 Зависимость структуры течения от величины зазора . 64
2.8 Влияние конвекции на форму поверхности эпитакси-ального слоя 68
2.9 Сравнение с результатами двумерных расчетов 69
2.10 Анализ влияния параметров сетки на результаты расчетов 71
2.11 Основные результаты 72
Математическое моделирование технологических режимов выращивания ЭС с предварительным подрастворением подложки 74
3.1 Описание процесса 74
3.2 Результаты расчетов 77
3.3 Учет конечной теплопроводности стенок ростовой камеры и жидкости 91
3.4 Основные результаты 95
Выводы
- Численный метод решения уравнений Навье-Стокса для несжимаемой жидкости
- Метод решения задачи концентрационной конвекции
- Влияние конвекции на форму поверхности эпитакси-ального слоя
- Учет конечной теплопроводности стенок ростовой камеры и жидкости
Введение к работе
Свободная конвекция, т.е. движение жидкости, возникающее в поле силы тяжести, вследствие пространственной неоднородности плотности, вызванной неоднородностью распределения температуры и (или) концентрации, - распространенное в природе явление. Оно наблюдается в системах любого масштаба: от микроскопического, в живых организмах, до масштабов океанов, земной атмосферы, Вселенной. Немалую роль конвекция играет и в разнообразных технических устройствах.
Методы численного моделирования давно и успешно используются для изучения конвекции. При этом проблемы, которые возникают в науке и технике, служат постоянным источником новых задач, требующих для своего решения разработки и использования специальных алгоритмов, учитывающих особенности изучаемых процессов. Одной из областей где, благодаря созданию эффективных вычислительных алгоритмов, математическое моделирование превратилось в важнейший инструмент исследований, является технология получения современных, в том числе полупроводниковых материалов.
В данной работе изучается концентрационная конвекция, как вид гидродинамической неустойчивости, применительно к проблеме получения полупроводниковых материалов методом жидкофазовой эпитаксии.
Метод жидкофазовой эпитаксии (ЖФЭ) - это способ выращивания монокристаллических слоев из раствора - расплава на подложке определённой кристаллографической ориентации [1]. В простейшем случае получения плёнок тройного полупроводникового соединения подложка состава AXoBi_XoC приводится в контакт с раствором компонентов А и В в расплаве С. Если при температуре контакта Т = То раствор - расплав является насыщенным, то жидкая и твёрдая фазы находятся в квазирав-
новесии. Для того чтобы вызвать процесс кристаллизации, температура системы расплав - подложка понижается с заданной скоростью. С понижением температуры падает растворимость компонентов А и В в С, раствор становится пересыщенным и из него на подложке кристаллизуется соединение AxBi_xC, х = х(Т), х(7о) = хо. В используемой модели предполагается, что в процессе роста эпитаксиального слоя на границе раздела фаз сохраняется квазиравновесие между раствором - расплавом и эпитаксиальным слоем, т.е. состав жидкой фазы на границе расплав -кристалл удовлетворяет фазовой диаграмме рассматриваемой системы. Вместе с тем основной объём расплава остаётся пересыщенным. Рост эпитаксиального слоя приводит к уменьшению концентрации растворенных компонентов в жидкой фазе вблизи растущего слоя. В результате формируется градиент плотности, который может быть как устойчивым, так и неустойчивым. Неустойчивый градиент плотности при превышении некоторого критического значения приводит к формированию конвективного движения. Перенос растворённых компонентов к фронту кристаллизации осуществляется механизмом конвекции и диффузии. Характер поступления растворённых компонентов к поверхности растущего слоя определяет скорость роста, равномерность её распределения вдоль поверхности эпитаксиального слоя (ЭС) и состав кристаллизующегося материала.
Хорошо известно, что конвективное движение в расплаве нарушает однородность распределения толщины пленки и ее состава по площади подложки [1]. С точки зрения качества получаемого материала, предпочтительными являются технологические режимы, обеспечивающие диффузионный режим роста или близкий к нему. Конкретные особенности массопереноса в растворе - расплаве зависят от условий технологического процесса: объёма расплава, скорости его охлаждения, положения подложек (вертикальное, горизонтальное, над расплавом или под ним), величины и направления температурных градиентов и т.д.
В данной работе рассматриваются системы в которых эпитаксиаль-ный рост осуществляется на подложку, расположенную горизонтально над расплавом или под ним. Изучаются реальные технологические режимы, в том числе такие, где процессу роста предшествует процесс рас-
творения подложки, в результате которого происходит удаление поверхностных дефектов и очищение подложки. Важнейшей задачей является определение критических значений параметров системы, при которых неподвижное состояние жидкости становится неустойчивым.
Численное моделирование процессов конвективной диффузии в растворе-расплаве при жидкофазовой эпитаксии ранее основывалось преимущественно на двухмерных моделях. Целью данной работы является изучение в трехмерной постановке гидродинамических процессов протекающих в расплаве при ЖФЭ, определение порога устойчивости системы, динамики развития неустойчивостей, анализ влияния возникающего конвективного движения на характеристики ЭС. Также исследуется зависимость наблюдаемых в расчетах течений и характеристик эпитаксиальных слоев от неоднородности распределения температур в ростовой камере. Проводится сравнение результатов, полученных в рамках двумерного и трехмерного приближений, определяется диапазон параметров, в котором существенно использование трехмерных моделей.
Процесс развития конвективной неустойчивости рассматривается в приближении Обербека-Вуссинеска [2, 3]. Соответствующая система уравнений получается из уравнений Навье-Стокса в предположении, что плотность жидкости не зависит от давления и является линейной функцией температуры и концентрации. Изменение плотности учитывается только в уравнении импульса в слагаемом, связанном с подъемной силой; градиенты скорости достаточно малы и тепловыделением за счет вязкой диссипации можно пренебречь. В результате процесс описывается системой уравнений Навье-Стокса для несжимаемой жидкости с подъемной силой в уравнении для импульса, и уравнениями конвективной теплопроводности и диффузии. Важную роль при моделировании процессов конвективного тепломассопереноса играет выбор алгоритма решения уравнений Навье-Стокса для несжимаемой жидкости.
О методах численного решения уравнений Наеье-Стокса для несжимаемой жидкости.
На протяжении многих лет усилия исследователей направлены на построение быстрых и надежных методов численного решения уравнений Навье-Стокса, однако до настоящего времени моделирование процессов на основе этих уравнений остается трудоемкой задачей, а свойства используемых алгоритмов изучены недостаточно. Особенно остро проблема выбора вычислительного алгоритма встает при моделировании длительных нестационарных процессов, когда метод решения должен обеспечивать получение физически корректных результатов на больших отрезках времени, желательно на сетках с небольшим числом узлов.
Характерной особенностью несжимаемой жидкости является бесконечная скорость переноса возмущений. Уравнение неразрывности не является эволюционными, в него входят лишь компоненты вектора скорости. В то же время уравнение, в котором в роли главной переменной выступает давление, - отсутствует. Это приводит к дополнительным трудностям при определении давления.
Для численного моделирования несжимаемых течений применяются два подхода. Первый из них состоит в том, что из исходной системы уравнений исключается давление и задача записывается в иных переменных: "функция тока-вихрь" [4,5], "вихрь-скорость" [4] или "функция тока-скорость" [6]. Второй подход основан на использовании естественных переменных скорость-давление [4,7-9].
Подход, основанный на использовании переменных "функция тока-вихрь", активно применяется для решения двумерных задач. Его достоинством является автоматическое выполнение условия несжимаемости, меньшее число неизвестных и отсутствие сложностей, связанных с определением давления. При этом, однако, возникают нелокальные граничные условия для вихря, строгая реализация которых требует совместного решения уравнений для вихря и функции тока или скорости [10-14]. Обобщение этого подхода на трехмерный случай сопряжено с немалыми трудностями и проигрывает алгоритмам в естественных переменных как по удобству использования, так и по числу уравнений [15]. Поэтому
в трехмерном случае основные усилия направлены на построение численных методов в переменных скорость-давление.
Дискретизация выбранной системы уравнений может осуществляться различными методами: спектральным, конечных разностей, конечных элементов и др. В данной работе используются конечно-разностные методы, обеспечивающие универсальность и простоту использования в областях сложной формы и при произвольных граничных условиях.
Основным требованием, предъявляемым к дискретной модели, является выполнение законов сохранения, справедливых в исходной физико-математической постановке задачи, так называемый принцип консервативности [16-18]. Уравнения Навье-Стокса для несжимаемой жидкости представляют собой записанные в дифференциальной форме законы сохранения массы (уравнение неразрывности) и импульса. Уравнение для кинетической энергии отсутствует, а баланс кинетической энергии является следствием выполнения законов сохранения массы и импульса. В разностном случае выполнение законов сохранения массы и импульса достигается аппроксимацией дивергентной формы исходных дифференциальных уравнений интегро-интерполяционным методом [17]. Однако численный метод должен сохранять не только массу и импульс, но и кинетическую энергию. Это особенно важно при моделировании существенно нестационарных задач на больших промежутках времени. В разностном случае из выполнения балансов массы и импульса, вообще говоря, не следует выполнения баланса кинетической энергии. Можно построить схемы, удовлетворяющие законам сохранения массы и импульса, но не сохраняющие кинетическую энергию [19,20].
Применительно к уравнениям Навье-Стокса для несжимаемой жидкости необходимость выполнения принципа консервативности впервые продемонстрировал A. Arakawa на примере алгоритмов в переменных "функция тока-вихрь" [21]. При интегрировании уравнений на больших временных отрезках нарушение в дискретной модели баланса кинетической энергии приводило к ее неограниченному росту. Аналогичное явление наблюдал Б.Л. Рождественский, используя гибрид метода конечных разностей и спектрального метода для моделирования турбулентности [22].
Если в дискретной модели закон сохранения энергии не выполнен, в системе наблюдается либо избыточная диссипация, либо нарастание кинетической энергии. В случае избыточной диссипации схема устойчива, но результаты, полученные с ее помощью при высоких значениях числа Рейнольдса, нефизичны: наблюдается сильное затухание энергии в областях высоких градиентов. Такими, например, являются схемы с аппроксимацией конвективных членов против потока первого порядка [19]. Нарастание кинетической энергии приводит к неустойчивости численной процедуры. Особенно опасны схемы, устойчивые при небольших значениях числа Рейнольдса и неустойчивые при больших. В расчетах с умеренным значением числа Рейнольдса они могут оставаться устойчивыми, но давать существенную погрешность за счет неконсервативности [19,23].
Для уравнений Навье-Стокса в переменных скорость-давления численные методы, сохраняющие кинетическую энергию, естественным образом строятся на разнесенных сетках [4,24]. В этих алгоритмах каждая компонента скорости и давление задаются на собственной сетке, конвективные члены аппроксимируются центральными разностями и схема имеет второй порядок точности по пространству. Существуют методы, использующие совмещенные сетки, когда все переменные рассматриваются в узлах одной и той же сетки. Это удобно с точки зрения программирования и при рассмотрении областей сложной формы в криволинейных координатах. Однако на совмещенных сетках приходится прибегать к дополнительным усилиям для сохранения кинетической энергии и для борьбы с сеточным осцилляциями давления. Они подавляются с помощью явного или неявного введения искусственной сжимаемости [25].
Сравнение алгоритмов на разнесенных и совмещенных сетках, выполненное в [26] для случая сжимаемых течений, показало, что на разнесенных сетках лучше выполняются законы сохранения. Другим преимуществом схем на разнесенных сетках является связь между скоростью и давлением, не дающая нефизических, сеточных осцилляции, характерных для схем на совмещенных сетках [20,27].
Схемы порядка точности выше второго и сохраняющие кинетическую энергию, можно найти в работах [19] и [20] для случая разнесенных и совмещенных сеток соответсвенно. Также были разработаны компактные
схемы высокого порядка точности [28]. S.K. Lele показал, что компактные схемы высокого порядка точности имеют меньший размер шаблона, и при определенном выборе параметров, лучше разрешают высокочастотные гармоники, чем стандартные конечно-разностные схемы того же порядка точности. Первоначально компактные схемы строились на основе формальной аппроксимации дифференциальных уравнений. В более поздних работах был применен интегро-интерполяционный метод [29]. Схемы высокого порядка точности, благодаря малой схемной вязкости, подвержены осцилляциям, приводящим к вычислительной неустойчивости [26,28,30] и, чем менее подробная сетка используется, тем выше вероятность их возникновения. Для борьбы с осцилляциями приходится вводить искусственную вязкость, фильтрующую высокочастотные гармоники, как, например, предложено в [28], или использовать регуляри-заторы высокого порядка точности [24,31-33].
Большинство методов решения уравнений Навье-Стокса для несжимаемой жидкости основано на дискретизации по времени типа предиктор-корректор [17,34,35]. Впервые алгоритмы такого типа были использованы в работах [36-38]. Их суть состоит в том, что на каждом временном слое сначала вычисляется предиктор скорости, вообще говоря, не удовлетворяющий условию несжимаемости. Затем из уравнения Пуассона определяется давление. При помощи найденного давления предиктор скорости корректируется таким образом, чтобы удовлетворить уравнению неразрывности. Это можно интерпретировать как проекцию предиктора скорости на множество соленоидальных функций, поэтому такие методы называют еще проекционными [37, 38]. Полученное после этапа коррекции поле скоростей считается найденным на текущем временном слое и осуществляется переход к следующему шагу по времени. В [36] использовалась разнесенная сетка, в [37,38] неразнесенная. Основные черты этих методов сохраняются до сих пор. Изложенные алгоритмы могут рассматриваться также как частный случай метода дробных шагов Яненко [34].
В проекционных методах на промежуточных этапах вычислений требуются краевые условия для предиктора скорости, в свою очередь определяющие краевые условия для давления. В 1972 С. Hirt и J. Cook [39]
предложили использовать метод [36] с уравнением Пуассона для определения поправки давления на следующем временном слое, вместо вычисления самого давления. Это позволяет достаточно просто ставить краевые условия для предиктора скорости со вторым порядком точности по времени [23]. Проблема выбора краевых условий обсуждается в [4,40], где показано, что краевые условия для давления, отсутствующие в дифференциальной постановке задачи, порождают ошибку, имеющую характер пограничного слоя, толщина которого пропорциональна корню из произведения вязкости на величину шага по времени [41].
В работах [41-45] показано, что порядок точности вычисления скорости на следующем временном слое совпадает с порядком аппроксимации метода расщепления. Точность определения давления при этом имеет лишь первый порядок по времени [46]. Использование явных формул, применявшихся в ранних методах в расчете предиктора скорости, приводит к серьезным ограничениям на величину шага по времени: 0.2bV2rRe < 1 и r/(Reh2) < 0.25 [4,47,48], избежать которых удается используя неявные и полунеявные схемы.
Существуют также алгоритмы, основанные на интегрировании уравнений Навье-Стокса по времени методом Рунге-Кутты [49] и Рунге-Кутты с совместным решением получающихся нелинейных уравнений итерациями по Ньютону [50]. Достоинством такого подхода является отсутствие необходимости ставить краевые условия для давления, что полностью согласуется с физической постановкой задачи. Однако решение нелинейных уравнений более трудоемко и сходимость таких методов часто вызывает затруднения. Попытка построить конечно-разностный метод, в которым система уравнений решается относительно вектора неизвестных, содержащего все компоненты скорости и давление, была предпринята в [48]. При решении двумерных задач, алгоритм оказался безусловно устойчивым, однако его обобщение на трехмерный случай связано со значительными трудностями, возникающими при решении соответствующей системы сеточных уравнений.
Для ускорения расчетов, ведутся разработки параллельных алгоритмов численного исследования процессов описываемых нестационарными уравнениями Навье-Стокса [51].
В данной работе для решения уравнений Навье-Стокса используется вычислительная процедура типа предиктор-корректор [7] второго порядка точности по пространству и первого по времени. Пространственная аппроксимация не вносит вклада в баланс кинетической энергии [24,52]. Метод использует центральные разности на разнесенной сетке. Давление относится к центрам ячеек, скорости - к центрам граней [53].
Неустойчивость Рэлея-Бенара
Начало систематическому исследованию конвективной неустойчивости положили экспериментальные работы Бенара. Он наблюдал возникновение движения в тонком, горизонтальном, подогреваемом снизу слое жидкости со свободной поверхностью [54]. Первые попытки теоретического исследования такой задачи были предприняты лордом Рэлеем, который провел линейный анализ задачи с двумя свободными границами [55]. Поэтому движение, возникающее в тонком горизонтальном слое жидкости под действием градиента температуры или концентрации, получило название конвекции Рэлея-Бенара. В настоящее время, этому вопросу посвящено огромное число работ, однако теория конвективной неустойчивости еще далека от своего завершения, ограничиваясь анализом частных задач и сопоставлением полученных результатов с экспериментальными данными [3,56-58].
Неустойчивость Рэлея-Бенара является одной из центральных проблем механики жидкости. Интерес к этому явлению объясняется тем, что в рамках простой физической постановки задачи удается получить широкий класс течений: от устойчивых и регулярных, до нерегулярных, быстро изменяющихся во времени. При этом в широком диапазоне параметров в толщине слоя сохраняется простая структура. Основные вопросы, возникающие при исследовании конвекции Рэлея-Бенара, это: изучение порога устойчивости рассматриваемой системы и планформы возникающего движения (проекции конвективной структуры на горизонтальную плоскость); анализ влияния структуры течения на характеристики тепломассопереноса; определение скорости развития конвективного движения из состояния покоя [3,58].
Экспериментальное исследование указанных явлений сопряжено с определенными трудностями. В первую очередь они связаны с большими временами эволюции конвективных структур (для воды от 10 минут до сотен часов, в зависимости от горизонтальных размеров области), высокой стоимостью экспериментальных установок, особенно для больших аспектных отношений, необходимых для полноценного анализа конвективных структур [59]. Шагом вперед в этой области было применение жидких газов, позволивших снизить характерные времена и расширить диапазон изменения значения числа Прандля [59]. Также представляет сложность само наблюдение за экспериментом и анализ структуры движения по полученным измерениям. Поэтому проведение численного моделирования и последующее сравнение полученных результатов с контрольными экспериментами является эффективным способом изучения конвективной неустойчивости.
Теоретические исследования устойчивости покоящейся жидкости проводятся, как правило, в стационарной постановке и в бесконечном горизонтальном слое. Температура верхней и нижней границ фиксирована, а внутри области изменяется линейно. Линейный анализ данной задачи позволяет найти критическое значение числа Рэлея. Критическому значению соответствует суперпозиция движений, имеющих форму параллельных валов определенной длины волны и произвольно ориентированных в горизонтальной плоскости. Амплитуда критического движения бесконечно мала. Для значений числа Рэлея выше критического, допустимым является целый диапазон волновых чисел [60-62]. В надкритической области любая допустимая мода нарастает неограниченно и найти нетривиальное стационарное решение в рамках линейного анализа не удается.
Найти и исследовать стационарные течения в надкритической области позволяет слабо-нелинейный анализ, основанный на представлении решения в виде ряда по степеням малого амплитудного параметра [3,63,64]. Этот параметр пропорционален корню квадратному из над-критичности: разности между значением числа Рэлея и его критической величиной. Сам по себе расчет стационарных надкритических движений не дает возможности ответить на вопрос, какое именно движение pea-
лизуется в действительности. В экспериментах могут наблюдаться лишь устойчивые, по меньшей мере к бесконечно-малым возмущениям, течения. Поэтому найденные решения исследуются на устойчивость при помощи линейного анализа. Наиболее полные результаты по исследованию устойчивости надкритических течений были получены А. Шлютером, Д. Лорцем и Ф. Буссе. Они показали, что в областях небольшой надкри-тичности устойчивыми относительно трехмерных возмущений являются только двумерные валы [64,65]. В дальнейшем были найдены диаграммы устойчивости различных типов течений, получившие название баллонов или сачка Буссе [62,65,66]. Однако механизм отбора планформ и волновых чисел до сих пор остается неопределенным [58].
Анализ влияния боковых стенок на структуру течения позволяет сузить спектр допустимых течений, сделав его дискретным. Но и в этом случае вопросы устойчивости и реализуемости течений остаются открытыми. При достаточно большом аспектном отношении влияние боковых стенок на критическое значение числа Рэлея мало, основное влияние стенок на структуру течения наблюдается в непосредственной близости от них [67-69]. Течение, как правило, имеет двумерную структуру, валы выстраиваются параллельно короткой стороне области или подходят к границам практически под прямым углом [68-70].
Техника проведения слабо-нелинейного анализа используется для изучения влияния температурной зависимости свойств жидкости, таких как вязкость или теплопроводность, на устойчивость надкритических течений [61]. В этом случае решение представляется в виде разложения в ряд по амплитуде и по малым параметрам, характеризующим зависимость свойств жидкости от температуры. Как показано в [61] для бесконечного числа Прандтля, в жидкости с переменной вязкостью возможно существование конечно-амплитудной подкритической конвекции, имеющей форму шестиугольных ячеек. В надкритической области, с ростом числа Рэлея шестиугольники становятся неустойчивыми относительно возмущений в форме валов. Направление движения жидкости в ячейках определяется характером зависимости параметров среды от температуры.
Более важным для рассматриваемого в данной работе класса задач
является случай, когда верхняя и нижняя границы области равномерно охлаждаются с постоянной скоростью. В результате внутри области формируется нелинейный профиль температуры. Эта задача эквивалентна стационарной задаче об устойчивости неподвижной жидкости в бесконечном горизонтальном слое с параболическим профилем температуры [71]. Нелинейный анализ, основанный на разложении в ряд по степеням малых параметров: амплитудного и кривизны профиля температуры, показал, что критическое значение числа Рэлея в данном случае ниже, чем при постоянном градиенте, и зависит от кривизны профиля [71]. Как и при переменной вязкости, здесь возможно возникновение конечно-амплитудной подкритической конвекции в форме шестиугольных ячеек, направление движения жидкости в которых определяется знаком кривизны профиля температуры. Эти задачи объединяет отсутствие симметрии верх-низ; для них характерны шестиугольные ячейки, в то время как для задач симметричных относительно средней плоскости, например при линейном профиле, типичны валы.
Как уже было сказано, практически все теоретические исследования устойчивости проводятся в стационарной постановке. Вопросы развития конвективного движения из начальных возмущений, скорость нарастания амплитуды различных форм течений теоретически почти не рассматриваются. Выбор начальных состояний, приводящих к течениям того или иного типа, определение характерных масштабов получающихся в результате течений - является самостоятельной, сложной и интересной задачей.
Для изучения реальных нестационарных задач со сложными краевыми условиями в областях произвольной формы приходится применять методы математического моделирования, основанные как на прямом численном решении уравнений Навье-Стокса, так и на решении приближенных уравнений, возникающих в ходе теоретических исследований.
Численное моделирование конвекции в трехмерной постановке началось достаточно давно. В [72] исследуется процесс формирования конвективных структур на основе двумерной модели трехмерной конвекции; описываются процессы возникновения и взаимодействия конвективных структур, а также их дефектов. В [73] для областей с небольшим ас-
пектным отношением (<3) в приближении бесконечного числа Прандтля изучается устойчивость двумерных валов относительно трехмерных возмущений. Численное исследование устойчивости жидкости с вязкостью, зависящей от температуры, для бесконечного значения числа Прандтля проведено в [74], где, в частности, было получено движение, охватывающее только часть слоя, с планформой в виде квадратных ячеек. Важность учета трехмерности показана в [15] на примере задаче о тепловой гравитационной конвекции в области с боковым подогревом, и в [75,76]. С увеличением значения числа Рэлея течение теряет двухмерную структуру и становится трехмерным. В [77] проводится численное моделирование подкритической конвекции в жидкости с вязкостью, зависящей от температуры. Для этого в шестиугольной планформе выбирается одна ячейка, ее внутренняя часть приближается цилиндром. На границе цилиндра ставятся периодические краевые условия и, используя предположение осесимметричности, задача сводится к двумерной. В работе получена зависимость критического числа Рэлея от перепада вязкости в области. В последнее время ведется интенсивное исследование конвекции Рэлея-Бенара в вязком сжимаемом газе [78-80]. Также, созданы пакеты программ, позволяющие осуществлять моделирование процессов тепломассообмена [81].
В данной работе изучается процесс потери конвективной устойчивости раствора-расплава при выращивании трехкомпонентных монокристаллических слоев методом ЖФЭ. Процесс возникновения и развития конвекции в системе во многом аналогичен случаю, рассмотренному в [71], но существенно нестационарен, т.к. граничные условия зависят от времени. На основе анализа результатов двумерных расчетов найдено значение критического числа Рэлея. В трехмерных расчетах зарегистрировано существование подкритической конвекции в форме шестиугольных ячеек; определено минимальное значение числа Рэлея, при котором наблюдается конвективное движение. Полученные результаты хорошо согласуются с теоретическими данными [71].
Учет теплопроводности степок.
Влияние конечной теплопроводности стенок на величину критического числа Рэлея было отмечено еще в 1926 году Н. Jeffreys [82]. В 1964 Е.М. Sparrow с соавторами выполнили анализ процесса развития неустойчивости в бесконечном горизонтальном слое, теплопередача на верхней и нижней границе которого определялась линейным законом Фурье, т.е. краевыми условиями третьего рода. Для нескольких случаев было определено критическое значение числа Рэлея в виде функции числа Био (характеризующего отношение теплопроводности внутри области и во внешней стреде) [83]. D.T. Hurle и Г.З. Гершуни с соавторами изучали устойчивость равновесия бесконечного горизонтального слоя жидкости, ограниченного сверху и снизу полубесконечными стенками, теплопроводность которых различна и отличается от теплопроводности жидкости [3,84]. Было показано, что чем ниже теплопроводность стенок, тем ниже порог устойчивости: он опускается с 1708 для идеально проводящих стенок до 720 для теплоизолированных. Волновое число при этом уменьшается с 3.11 до нуля.
Исследованию влияния боковых стенок на порог устойчивости посвящены работы [67,85-88]. Основные результаты этих исследований состоят в следующем: теплопроводность боковых стенок влияет на значение критического числа Рэлея только при аспектном отношении меньше двух, при этом значение критического числа Рэлея для области с непроводящими стенками значительно ниже, чем в случае идеально теплопро-водящих боковых границ. Очевидно, что первый результат обусловлен вязким взаимодействием жидкости с боковыми стенками, а второй - тем что температурные флуктуации вблизи идеально теплопроводной стенки затухают быстрее, чем возле теплоизолированной. Эксперименты, проведенные К.Stork и и.МШІег'ом для области с различной теплопроводностью верхней и нижней стенок и с плохо проводящими боковыми границами подтвердили эти результаты [68].
Работы по изучению влияния конечной теплопроводности боковых стенок на значение критического числа Рэлея и реализующиеся формы течений продолжаются и в настоящее время, однако вопрос в целом все
еще остается открытым. Из более поздних работ отметим [89], в которой рассматривается бесконечный горизонтальный слой жидкости, ограниченный стенками конечной толщины и теплопроводности. Расчет линеаризованной задачи показал, что при отношении толщины горизонтальных границ к толщине жидкой фазы больше единицы, результаты близки к полученным Г.З. Гершуни для полубесконечных стенок [3]. Если же это отношение порядка единицы, то для стенок с теплопроводностью на порядок выше чем в жидкости, разница между двумя результатами достигает 20%. Для волнового числа это влияние более существенно. При отношения толщин порядка единицы, разница между результатами [89] и [3] в зависимости от теплопроводностей стенок и жидкости может достигать 20%, если стенки в 10 раз тоньше слоя жидкости - то разница может возрастать до 100%. В [89] также отмечалось, что даже небольшие перепады температуры у границы, являются причиной возникновения пристеночного вихря.
Таким образом, учет конечной теплопроводности и толщины стенок может оказать существенное влияние на наблюдаемые картины течений в задаче Рэлея-Бенара. В данной работе проводится численное исследование процесса роста кристалла методом ЖФЭ с учетом конечной теплопроводности стенок ростовой камеры. Результаты расчетов сравниваются с полученными для случая идеально-проводящих стенок.
Содерэюание диссертации по главам
Во введении содержится общая характеристика рассматриваемого класса задач, обзор литературы, посвященной методам численного исследования задачи, а также известные результаты экспериментальных, аналитических и численных исследований, необходимые для лучшего понимания свойств описываемого класса задач и целей работы.
Первая глава посвящена построению и анализу алгоритмов численного исследования процессов конвективного тепломассообмена. Основу математической модели составляют уравнения Навье-Стокса в приближении Буссинеска, уравнения конвективной теплопроводности и диффузии. Задача рассматривается в трехмерной постановке в естествен-
ных переменных. При построении алгоритма используется процедура расщепления по физическим процессам. На каждом шаге по времени, сначала решаются уравнения Навье-Стокса, а затем, по найденному полю скоростей определяется распределение температуры и концентрации. Дискретизация по времени имеет первый порядок точности.
Для решения уравнений Навье-Стокса используется известный полунеявный метод предиктор-корректор на разнесенных сетках. Его суть состоит в замене уравнения несжимаемости уравнением Пуассона для давления. Сначала из уравнений движения вычисляется предиктор скорости. Поле скоростей на этом этапе, вообще говоря, не удовлетворяет условию несжимаемости. Потом находится поправка давления, и с ее помощью осуществляется коррекция поля скоростей для обеспечения несжимаемости. Пространственная аппроксимация уравнений выбирается исходя из требований консервативности и имеет второй порядок внутри области и первый на границе.
В работе проводится анализ различных способов решения разностного уравнения Пуассона для давления. В первом из них используются краевые условия Неймана, т.е. задача для давления является вырожденной. Ее можно решать методом сопряженных градиентов, находя решение наилучшего приближения. В других способах, значение давления фиксируется в одной точке и получается невырожденная задача. Проведенный теоретический анализ показал, что условие несжимаемости нарушается в той ячейке сетки, где фиксируется давления. Величина отклонения от несжимаемости равна суммарной невязке решения уравнения, умноженной на шаг по времени и деленной на объем (или площадь) ячейки. Таким образом, использование не неймановских краевых условий может приводить к серьезным нефизическим ограничениям на шаг по времени и пространству. Этот вывод подтверждается результатами численных расчетов с двойной и одинарной точностью в двумерном и трехмерном случаях. Кроме того, число итераций в методе сопряженных градиентов при использовании краевых условий Неймана, значительно ниже, чем для других рассмотренных случаев.
Тестирование алгоритма решения уравнений Навье-Стокса проводилось в двумерном и трехмерном приближениях на примере задачи о
движении жидкости в каверне с движущейся крышкой, и на задаче о естественной тепловой конвекции. Тестирование показало надежность и устойчивость используемого алгоритма в широком диапазоне параметров.
Во второй главе приводится математическая постановка задачи о выращивании тройных монокристаллических соединений методом жидко-фазовой эпитаксии. Основной особенностью задачи является наличие нелинейных краевых условий на границе расплав-кристалл, связывающих между собой значение концентрации растворенных компонентов в жидкой и твердой фазах. В работах [90-92] показано, что численная реализация подобных условий требует совместного решения уравнений для концентрации на границе расплав-кристалл. Дискретизация системы уравнений по времени и метод решения уравнений Навье-Стокса были приведены в первой главе. В данной главе осуществляется пространственная аппроксимация уравнений конвективной диффузии с учетом требований консервативности по концентрациям. Для решения нелинейных на границе раздела фаз уравнений используется метод Ньютона.
С точки зрения гидродинамики, рассматриваемая задача является задачей о развитии конвективной неустойчивости в неподвижной жидкости с нестационарными и нелинейными краевыми условиями для концентрации. В работе строится последовательность приближений, позволяющих использовать известные теоретические и экспериментальные данные по гидродинамической неустойчивости для анализа результатов расчетов.
Далее, проводится численное моделирование процесса роста ЭС в трехмерной постановке для значений числа Рэлея в диапазоне 1.1-103 — 4.5-Ю5. Изменение значения числа Рэлея осуществлялось за счет выбора толщины слоя жидкости. При моделировании обнаружена подкрити-ческая конвекция, которая является существенно трехмерной, и имеет форму шестиугольных ячеек. Форма подкритического движения и минимальное значение числа Рэлея, при котором оно наблюдается, полностью согласуются с результатами конечно-амплитудной теории устойчивости [3,58,71]. В надкритической области, для значений числа Рэлея в диапазоне 1.6.. < Ra < 5Racn Rao- = 1545 зарегистрирована устой-
чивая ячейковая конвекция, при 5Racr < Ra < 2(Шасг спицевидная, и для Ra ^ 2QRacr хаотическая.
Кроме того, изучено влияние толщины слоя жидкости на наблюдаемую в результате развития неустойчивости структуру течения и на форму поверхности выросшего кристалла. Неоднородность поверхности выросшего слоя определяется структурой течения, его устойчивостью и амплитудой, поэтому важное значение имеет множество допустимых течений, которое в трехмерном случае значительно шире, чем в двухмерном. Сопоставление характеристик эпитаксиальных слоев, полученных в трехмерном случае, с данными двумерных расчетов показало, что:
существенно трехмерная конечно-амплитудная подкритическая конвекция приводит к значительной неоднородности слоя там, где в двумерной модели наблюдается диффузионный режим;
в областях слабой надкритичности, в трехмерном случае наблюдается более интенсивная конвекция, приводящая к большей неоднородности распределения толщины эпитаксиального слоя по площади подложки, чем двухмерном;
в областях большой надкритичности, скорость перестройки трехмерных течений выше, а неоднородность поверхности выросшего слоя - меньше, чем в двумерном.
В конце главы проводится анализ влияния величины шагов сетки по времени и пространству на результаты моделирования. В описанных выше расчетах потеря устойчивости происходит под действием возмущений, порождаемых самим численным алгоритмом. Они формируют структуру конвективного движения на начальной стадии процесса, от которой зависит дальнейшая эволюцию течения. Измельчение сетки влияет на возмущения и, тем самым, на всю динамику процесса. Поэтому на один и тот же момент времени решения, полученные на различных сетках, могут заметно отличаться. Для сравнения таких решений необходимо использовать интегральные характеристики: среднюю кинетическую энергию, характерный размер и форму конвективных структур. Анализ показал, что указанные характеристики решений слабо зависят от параметров используемых сеток.
В третьей главе изучаются режимы эпитаксиального выращивания с предварительным растворением подложки. Подложка, расположенная в верхней части ростовой камеры, приводится в контакт с ненасыщенным раствором и температура всей системы понижается. В результате начинается процесс растворения и в жидкости формируется неустойчивый градиент плотности. Через некоторое время происходит насыщение жидкой фазы материалом подложки и начинается процесс эпитаксии, в ходе которого градиент плотности становится устойчивым.
Анализ влияния конвективного движения, возникающего в ходе растворения, на свойства эпитаксиального слоя начинается с рассмотрения процесса в изотермических условиях. В этом случае неустойчивый градиент плотности сохраняется дольше, чем в охлаждаемой системе, что позволяет течению набрать большую кинетическую энергию и сформировать четкую структуру. Чем выше интенсивность конвективного движения и устойчивее его структура, тем выше неоднородность фронта травления. Таким образом, мы получаем задачу, для которой время существования и развития конвекции является важнейшей характеристикой.
Существенным фактором, определяющим время развития конвекции и неоднородность получаемого слоя, является величина и форма возмущений, вызывающих потерю устойчивости. Возмущения, вносимые численным методом, трудно контролировать, поэтому была проведена серия расчетов с заданным начальным возмущением в виде двумерных валов или трехмерных квадратных ячеек. Результаты расчетов для различных значений толщины жидкой фазы и степени недосыщения расплава показали, что скорость роста амплитуды двумерных валов выше, чем у трехмерных ячеек, и они дают большую неоднородность фронта травления. Максимальная неоднородность, тем не менее, достигается при сочетании в сформировавшемся течении валов и ячеек. Двумерные валы быстро набирают энергию и передают ее трехмерным структурам, которые обеспечивают более высокую неоднородность распределения потока растворенных компонентов по площади подложки.
Роль начального возмущения особенно велика в задачах с небольшой надкритичностью. Без начального возмущения слой получается более
ровным, так как конвекция не успевает набрать достаточную энергию. С ростом надкритичности начальное возмущение перестает играть значимую роль в формировании характеристик слоя, хотя скорость нарастания кинетической энергии по прежнему максимальна для двумерных валов. Поэтому моделирование процесса в областях большой надкритичности можно проводить с однородными начальными данными. Результаты расчетов на измельченной сетке показывают, что начальное возмущение оказывает большее влияние на неоднородность фронта травления, чем сеточные параметры.
В двумерном случае имеет место аналогичная картина: конвекция развивается быстрее при внесении начального возмущения, однако когда надкритичность велика, это не влияет на неоднородность слоя. Двумерные и трехмерные результаты достаточно близки, поскольку в трехмерном случае наибольшей скоростью роста среди рассматриваемых возмущений обладают именно двумерные структуры. Таким образом, с прикладной точки зрения при изучении процесса растворения надежные результаты обеспечиваются в рамках двумерного подхода.
Численное моделирование процесса роста в неизотермических условиях подтверждает очевидный факт: время существования конвекции и неоднородность фронта травления уменьшаются при увеличении скорости охлаждения системы. Выводы, сделанные относительно влияния начальных возмущений на динамику процесса, справедливы и в этом случае.
Далее, постановка задачи усложняется введением конечной теплопроводности и конечной толщины стенок ростовой камеры. Математическая модель дополняется уравнением теплопроводности, которое решается в расплаве и внутри стенок камеры. В расчетах используется разностная схема с консервативной аппроксимацией конвективных членов. В технологических режимах выращивания ЭС, скорость охлаждения мала, вместе с тем теплопроводность расплава, из которого осуществляется рост, велика, поэтому распределение температуры в жидкости практически однородно по пространству и задача близка к ранее рассмотренному случаю идеально проводящих среды. Ранее предполагалось, что неоднородностью температуры внутри ростовой камеры можно пренебречь, по-
скольку при рассматриваемых технологических режимах значение температурного числа Рэлея на три порядка меньше концентрационного, и конвективное движение, вызываемое неоднородностью температуры, не влияет на динамику роста и растворения. Однако результаты расчетов показывают, что неоднородность фронта травления выше при конечной теплопроводности. Объясняется это тем, что горизонтальная неоднородность температуры играет роль постоянно действующего возмущения фиксированной формы, в отличие от вносимого только в начальный момент времени возмущения в рассматривавшейся ранее задаче с идеально проводящей средой. Такое возмущение увеличивает скорость развития конвективного движения и приводит к росту неоднородности слоя.
Основные результаты работы:
Разработан разностный алгоритм и комплекс программ для численного исследования трехмерных нестационарных процессов конвективного тепломассопереноса применительно к проблеме выращивания трехкомпонентных эпитаксиальных слоев из жидкой фазы.
Впервые проведено численное исследование процесса выращивания тройных полупроводниковых соединений из жидкой фазы в трехмерном приближении с учетом реальной фазовой диаграммы системы. В расчетах обнаружена подкритическая конвекция в форме шестиугольных ячеек. В надкритической области, в зависимости от значения числа Рэлея, зарегистрирована ячейковая, спицевид-ная и хаотическая конвекция. Результаты расчетов согласуются с данными линейного и слабо нелинейного анализа устойчивости.
В диапазоне изменения параметров, представляющем практический интерес, найдены технологические режимы, при которых мас-соперенос в жидкой фазе осуществляется преимущественно механизмом диффузии. В системах с предварительным растворением подложки проанализировано влияние динамики развития возмущений полей температуры и концентрации на характеристики ЭС. Показано, что скорость роста амплитуды возмущений в форме дву-
мерных валов выше, чем у трехмерных ячеек, и они приводят к большей неоднородности поверхности слоя. Перепад температуры порядка 0.1-0.3 градуса при абсолютной величине Т ~ 700ІІГ, порождаемый конечной теплопроводностью веществ, приводит к увеличению неоднородности толщины выросшего слоя по сравнению с идеально теплопроводящей средой в 1.5-2 раза.
Автор выражает искреннюю благодарность:
своему научному руководителю д.ф.-м.н. Мажоровой О.С. и член, корреспонденту РАН, д.ф.-м.н. Попову Ю.П. за постоянное внимание и помощь в работе, участие в многочисленных обсуждениях; сотрудникам 27-й лаборатории института ГИРЕДМЕТ за постановку задачи; всем сотрудникам 11-го отдела ИПМ им М.В. Келдыша РАН за сотрудничество и ценные советы во время подготовки диссертационной работы.
Численный метод решения уравнений Навье-Стокса для несжимаемой жидкости
Опыт численного решения уравнений Навье-Стокса в естественных переменных показывает, что способ нахождения давления на новом временном слое в значительной степени определяет точность и устойчивость всей вычислительной процедуры. Построенная в предыдущем параграфе схема предполагает решение задачи Неймана для давления, оператор которой является вырожденным. Его ядро состоит из констант, и все множество решений для конкретной правой части может быть получено сложением частного решения с произвольной постоянной. Условием разрешимости задачи Неймана является ортогональность правой части уравнения (в рамках этого параграфа обозначим ее F) ядру, т.е. (F, 1)2/1 = 0 [95, стр. 479]. В нашем случае F = VhV/т и это условие выполнено. Однако при численной реализации, из-за погрешностей вычисления правой части, это условие выполнено уже не будет.
Для однозначной разрешимости (1.3.9)-(1.3.10) потребуем, чтобы в некоторой точке А Є 0,р и находящейся на границе области, давление равнялось нулю, т.е. РА = &РА — 0. Реализовать это можно несколькими способами: N. Методом сопряженных градиентов найдем др, минимизирующее функционал Ддр — F\\2h- др существует для любой правой части F и является решением задачи Неймана А др = F — (F, 1)2 . Искомое решение получается следующим образом: др = др — дрА. NP. Потребуем, чтобы в точке А вместо (V dp, га) = 0 было выполнено 5рА = 0. При таких краевых условиях задача однозначно разрешима для любой правой части. Но выполнено ли при этом условие несжимаемости? Проверим это, используя для наглядности рис. 1.1. Выберем точку А в центре приграничной ячейки, содержащей А, и найдем значение дивергенции скорости в этой ячейке сетки. Пусть F(i,j,k) - погрешность вычисления правой части. Проинтегрируем уравнение (1.3.9) по всей области: (Aph6p, lb = (F + F, Цгл = (F, 1)2 (Aph6p,i)2h= Е Д№А = - Е (vs ,7g(ds,n)M Учитывая что во всех граничных узлах, кроме А, задано условие (V dp, п) = 0, получаем: (Vphdp,n)\(dS,n)\A = -(F,l)2h Граничное условие для нормальной компоненты скорости в точке А - У!1 = п,А запишем в виде, аналогичном (1.3.11): К/1 = УпЛ - r(Vphdp, п)А + r(Vphdp, п)А (1.4.1) Используя уравнения (1.3.11) и (1.4.1), преобразуем: Vfc17+1 = VhVA - rAph5pA - r(VphSp,n)A/(hA,nA) = = {VphSp,n)A/(hA,nA) \VhVA\ = \r(F, l)2h/dnA\ (1.4.2)
Пусть решение задачи для давления (1.3.9)-(1.3.10) найдено и невязка ty=Apl8p—VhY/т не превосходит по модулю некоторой заданной величины є/r. Естественно потребовать, чтобы значение IV V 1] из (1.3.11)-(1.3.12) было не больше е. При использовании краевых условий Неймана (N) во всех точках области V/jVm+1=T\I/, и это требование выполняется. Для остальных двух способов (NP), (NPP) равенство УьУт+1=тФ справедливо всюду, кроме точки А, в которой.
Найдем условия, при которых величина (1.4.5) остается достаточно малой. Значение (F, 1)2/1, вообще говоря, зависит от времени, определяется точностью выполнения арифметических операций и порядком самой правой части F. В расчетах задачи о движении жидкости в каверне при Re « 10 -f-103 с одинарной точностью (F, 1)2 Ю-8, для двойной точности (F, l)2/i 10"14. Пусть порядок величины (F, 1)2/1 сохраняется при измельчении шага сетки по пространству и по времени, и выясним, при каких условиях уравнение несжимаемости будет выполнено с точностью є = Ю-5.
В двумерном случае, при одинарной точности, оценка приобретает вид: т10_8//і2 10"5, или т 103/i2. Аналогично получаются условия для двойной точности и в трехмерном случае. Порождаемые этими условиями ограничения на величину г, вычисленные для различной величины шага сетки по пространству, приведены в следующих таблицах: Таблица 1.1: Двумерный случай шаг по пространству К = Ю"1 К = Ю"2 К = 10"3 одинарная точность т 10 г Ю-1 г 10"3 двойная точность г 107 г 105 т 103 Таблица 1.2: Трехмерный случай шаг по пространству К = Ю-1 /г = 10 2 /і = 10"3 одинарная точность т 1 т 10"3 г 10"6 двойная точность г 106 т 103 г 1
Данные, приведенные в таблицах показывают, что в двумерном случае на сетках со сравнительно небольшим числом узлов использование краевых условий (NP), (NPP) при одинарной точности еще не влечет за собой серьезных ограничений на величину шага по времени. Вместе с тем в трехмерном случае применение краевых условий (NP), (NPP) оправдано только при использовании двойной точности.
Однако оценка (1.4.2) лишь качественно описывает зависимость точности выполнения условия несжимаемости от параметров расчета. Реальные ограничения на величину шага по времени, обусловленные требованием выполнения условия несжимаемости с необходимой точностью, еще более жесткие. Так расчеты с двойной точностью двумерной задачи о течении жидкости в каверне на сетке с шагами h\ = fi2 = h = 0.05, г =1 и при є = 10 5, показывают, что для краевых условий Неймана условие несжимаемости diVhV = 0 выполняется с точностью 5 10 6, и с точностью до 4 Ю-4, 1 10 4 для краевых условий (NP), (NPP) соответственно. Для (NP), (NPP) максимальное значение дивергенции достигается в точке А [53].
В данной серии расчетов значение коэффициента а в предикторе давления выбиралось равным 1, что обеспечивает аппроксимацию урав -39 нения (1.1.1) уравнением для предиктора скорости (1.3.2). Для решения систем линейных алгебраических уравнений используется итерационный метод Ortomin [96-98]. Он относится к классу методов сопряженных градиентов с матрицей предобуславливания, построенной с помощью неполного LU разложения исходной матрицы. Итерации продолжались до тех пор, пока величина невязки не становилась меньше Ю-5. Условием достижения стационара является выполнение неравенства S( m + 20) - S(tm)\\c 5 10"5, где S=(V,P).
Влияние выбора граничных условий для давления на сходимость итерационного процесса. На рис. 1.2 для трех рассматриваемых способов задания граничных условий для давления изображена зависимость числа итераций в методе Ortomin от номера шага по времени при решении задачи о движении жидкости в каверне с Re = 1000. Расчеты проводились с одинарной точностью на сетке 26 х 14 с шагом по времени т = 1.
Метод решения задачи концентрационной конвекции
1. Проведен анализ влияния способа задания граничных условий в уравнении Пуассона для давления на свойства вычислительной процедуры. Показано, на одинаковых сетках в расчетах с граничными условиями Неймана точность выполнения условия несжимаемости на порядок выше, а численная реализация методом Ortomin эффективнее по числу итераций в 2-3 раза, чем при использовании других рассмотренных типов граничных условий.
2. Проведено тестирование алгоритма (1.3.7)-(1.3.12) на решении модельной задачи течения жидкости в квадратной каверне с движущейся крышкой для Re = 10 4- 1000. Сравнение результатов полученных на сетках 20x20, 50x50 с прецизионными расчетами на сетке 256x256 [7], показывало отклонение не более 10%.
Аналогичное сравнение для случая тепловой конвекции с неустойчивым градиентом плотности при Ra 20000 [11], дало отклонение не более 7%. Возникающее при этом ограничение на шаг сетки по времени для стационарной задачи свободной конвекции слабо зависит от величины шага сетки по пространству.
3. Для задачи о движении жидкости в каверне с движущейся крыш кой проведено сравнение течений в трехмерной и двумерной кавер нах, и найдена область вида: х% Є [ує, L 2 — ує] в которой отличие между ними не превышает заданной величины є, т.е. в этой области вместо решения трехмерной задачи можно использовать решение двумерной задачи, и отклонение будет не более чем на е.
Рассмотрим процесс выращивания тройных эпитаксиальных соединений AxBi_xC из раствора компонентов А и В в расплаве компонента С. Пусть ростовая камера, имеющая форму прямоугольного параллелепипеда (рис. 2.1), заполнена раствором-расплавом. На дно ее помещена подложка АХоВі_ХоС. Вся система находится в поле силы тяжести. В начальный момент времени расплав и подложка имеют одинаковую температуру (То) и находятся в квазиравновесии, т.е. концентрация растворенных компонентов и состав твердой фазы удовлетворяют фазовой диаграмме системы. Затем температура всей системы понижается по заданному закону. С понижением температуры падает растворимость компонентов А и В в расплаве С, раствор становится пересыщенным и начинается процесс его кристаллизации на подложке.
При построении математической модели процесса будем предполагать, что эпитаксиальный рост происходит в квазиравновесном режиме, т.е. концентрация растворенных компонентов на фронте кристаллизации и состав растущего слоя удовлетворяют фазовой диаграмме системы. Кинетические явления на границе раздела фаз не учитываются. В частности, это означает, что скорость роста эпитаксиального слоя определяется скоростью поступления растворенных веществ к фронту кристаллизации и не лимитируется процессами, протекающими на межфазовой границе. Перенос растворенных в расплаве компонентов осуществляется механиз -49 Рис. 2.1: Схема ростовой камерымами диффузии и естественной конвекции.
Возможность возникновения в расплаве конвекции связана с зависимостью плотности жидкой фазы от ее состава, температуры и давления. Однако зависимость плотности расплава от температуры в данном случае не существенна. В рассматриваемых системах ( см. таблицу А.1) коэффициент температуропроводности Л на 3-4 порядка больше коэффициентов диффузии . Характерный масштаб времени, связанный с процессом теплопроводности t\ = Щ/Х на порядок меньше, чем характерный масштаб времени вязкой диссипации tv = НЦи и на три порядка меньше характерных времен диффузии д = H /Di, где Н - характерная высота области. В результате пространственная неоднородность температуры, обусловленная охлаждением стенок ростовой камеры мала и не оказывает заметного влияния на величину подъемной силы. Таким образом, при "невысоких" скоростях охлаждения температуру расплава можно считать однородной по пространству и изменяющейся во времени по заданному закону Т = T(t).
Таким образом, основной причиной возникновения конвекции является зависимость плотности от состава. Рост эпитаксиального слоя при -50 водит к обеднению жидкой фазы вблизи растущего слоя по компонентам А и В. В результате формируется вертикальный градиент плотности, который является устойчивым, если плотность жидкости вблизи подложки больше, чем в верхней части расплава, и неустойчивым в противном случае.
В данной работе рассматривается процесс эпитаксиального роста из расплава, плотность которого падает с уменьшением концентрации растворенных компонентов. Так как рост осуществляется на подложку, расположенную на дне ростовой камеры, то формирующийся в процессе эпитаксии градиент плотности является неустойчивым и приводит к возникновению конвективного движения. Жидкость предполагается несжимаемой и гидродинамические процессы рассматриваются в приближении Буссинеска [2,3].
Влияние конвекции на форму поверхности эпитакси-ального слоя
Описанный процесс иллюстрируют рисунки Б.10 и Б.11. На рис. Б. 10 на последовательные моменты времени изображены графики зависимости среднего значения концентрации Cav в горизонтальных плоскостях от вертикальной координаты плоскости z. Из рисунка видно, что при t t\=bmin конвективное движение не оказывает существенного влияния на распределение концентрации. В момент t\ - формируется профиль концентрации, отчетливо свидетельствующий о наличии развитого конвективного движения в расплаве. На интервале времени ti t t2=10.5min происходит выравнивание распределения концентрации по высоте. При t t2 наблюдается медленное увеличение кривизны профиля.
В качестве характеристики неоднородности распределения концентрации в плоскостях, параллельных подложке, введем величину D(z) = C7(nr», yj, z) — Cav(z)\\Lhtx у График этой функции на различные моменты времени приведен на рис. Б. 11. Конвективное движение возникает в момент времени to 3mm, в этот же момент появляется пик на гра -66 фике функции D(z), однако его значение достаточно мало и график соответствующей кривой на рисунке не приводится. На момент времени ti=5min неоднородность становится максимальной, затем резко падает и продолжает уменьшаться вплоть до t2=10.5min. Затем D{z) монотонно возрастает во всех горизонтальных сечениях на протяжении всего времени расчета: t 80mm(10.6 ).
Возникающая на конечной стадии расчета картина течения приведена на рис. Б.12. Сравнивая рисунки Б.7 и Б.12 видно, что постепенная перестройка течения приводит к увеличению размеров конвективных ячеек. К моменту времени t = 80mm этот процесс еще не завершен. Наблюдаемые ячейки являются ячейками 1-типа.
При изменении толщины жидкой фазы в диапазоне от 1.5 -т- 2, как видно из рис. Б.8, графики зависимости значения числа Рэлея от времени ведут себя практически одинаково. Типичным является наличие резкого максимума, за которым происходит быстрое падение значения числа Рэлея с последующим медленным монотонным ростом. Аналогичная зависимость наблюдается и для кинетической энергии (рис. Б.9). Экстремальные значения Ekin(t) достигаются на более поздние, по сравнению с Ra(t), моменты времени, что обусловлено инерцией движения жидкости.
Течение имеет ячейковую структуру с характерным размером, зависящим от толщины жидкой фазы. На рассматриваемых временах процесс развития конвективного движения представляет собой последовательное "исправление" дефектов, существующих в структуре течения, за счет объединения ячеек в более крупные. На рисунках Б.13 и Б.14 показан процесс укрупнения ячеек для зазора Н = 2.0.
При дальнейшем увеличении величины зазора в диапазоне Н = 2 -f- 3 эволюция течения происходит по прежнему сценарию, однако графики Яа(), Ekin(t) приобретают колебательный характер (рис. Б.15, Б.16), обусловленный быстро протекающими процессами перестройки течения.
В случае Н = 2.5 течение уже не имеет четко выраженной ячеистой структуры (рис. Б.17, =5?шп(0.24 )). Оно похоже на сложную неустойчивую систему валов. С течением времени характерный размер возникающих структур увеличивается (рис. Б.18, і=АЬтгп{2.\6ів))- В экспе -67 риментальной работе [107] наблюдался похожий характер течения, при котором постоянно взаимодействуют, объединяясь и распадаясь, фрагменты валов. В указанной работе, такой характер течения связывался со спицевидной конвекцией, которая является развитой формой узелковой неустойчивости [58, с. 84-85]. В расчетах, в процессе развития течения, участки валов стремятся расположиться перпендикулярно границам области, что согласуется с экспериментальными и теоретическими данными приводимыми в параграфе 2.5. При больших значениях числа Rа спицевидная структура течения флуктуирует во времени и имеет хаотический характер [58, с. 134].
В расчетах хаотическое движение было получено при Н = 3, когда на развитой стадии движения число Рэлея становится равным Ra « 2ARacr. Вместе с тем наблюдавшаяся тенденция к увеличению характерного размера конвективных структур, связанная с эволюцией течения во времени и увеличением толщины жидкой фазы, нарушается. На ранних стадиях процесса характерный размер ячеек примерно совпадает с наблюдаемым при Н = 2 (рис. Б. 19, t= 3min(0.lti))), а затем уменьшается и начинается хаотическая конвекция (рис. Б.20, =45гат(1.5о)).
Расчеты проведенные для случая прямоугольной подложки со сторонами Li=25, 1/2=50, показали, что в рассматриваемом диапазоне значений величины зазора (Н = 1.5- 3) изменение горизонтальных размеров области практически не влияет на структуру течения на его развитой стадии. Характер зависимости числа Рэлея и кинетической энергии от времени сохраняется. Близки также и времена достижения функциями Ra(t) и Ekinif) своих экстремальных значений (рис. Б.21). Однако в случае прямоугольной подложки возникающие на начальной стадии расчета возмущения не всегда однородно распределены по пространству. В результате переход к конвекции осуществляется не одновременно по всей области, и происходит небольшое размывание экстремумов функций во времени, более заметное на больших зазорах.
В качестве примера приведем результаты, полученные для прямо -68 угольной подложки при Н = 1.5. Как и в случае квадратной подложки, течение зарождается на момент времени to — Зтт(0.4 ). Его начальная структура имеет форму ячеек (рис. Б.22). На развитой стадии конвекции, в момент максимума кинетической энергии i=5mm(0.67c ), структура течения имеет вид ячеек, большинство которых представляет собой шестиугольники неправильной формы. В процессе эволюции, происходит формирование ячеек правильной шестиугольной формы, что иллюстрирует рис. Б.23, соответствующий минимуму значения кинетической энергии при 2=10шт(1.3 )), и рис. Б.24 для =43тт(5.7о). Процесс эволюции к указанному моменту времени еще не завершен. Однако тенденция к установлению движения в виде правильных шестиугольников четко прослеживается, что согласуется с теоретическими результатами [71].
Представление о характере движения, наблюдавшегося в расчете при зазоре Н = 2.5, дает рис. Б.25, =5mm(0.24c ). Сравнение его с рис. Б. 17, где приведено распределение концентрации на верхней границе области на тот же момент времени для случая квадратной подложки, показывает, что течения имеют сходную структуру: представляет собой сочетание ячеек неправильной формы и фрагментов прямых и изогнутых валов.
Учет конечной теплопроводности стенок ростовой камеры и жидкости
Анализ безразмерных параметров задачи, проведенный в начале второй главы (стр. 49), показал, что значение теплового числа Рэлея много меньше диффузионного и меньше критического. Таким образом, неоднородность температуры, возникающая в расплаве, достаточно мала и не приводит к возникновению тепловой конвекции. Вместе с тем неоднородное поле температуры в ростовой камере приводит к неоднородности поля концентрации в горизонтальных плоскостях. В результате формируется возмущение специфической формы и жидкость выходит из состояния равновесия (нарушается необходимое условие равновесного состояния жидкости). Далее, в зависимости от величины надкритичности, это возмущение может нарастать. Поскольку характеристики конвективного движения существенно зависят от формы вносимого возмущения, рассмотрим задачу с конечной теплопроводностью стенок ростовой камеры и жидкости.
На границе жидкость-контейнер поток тепла остается непрерывен. В правой части уравнений Навье-Стокса теперь учитывается подъемная сила, связанная с неравномерным нагревом. Граничные условия для скорости и концентраций остаются прежние. На внешней границе контейнера Т = TQ — at, х Є dlext Дискретизация системы (3.3.1)-(3.3.4) проводится аналогично предыдущим случаям. Остановимся подробнее на построении разностного аналога уравнения (3.3.4).
Такая аппроксимация обеспечивает выполнение баланса тепловой энергии, а разностные аппроксимации операторов конвективных членов и теплопроводности, наследуют свойства дифференциальных аналогов, аналогичные (1.3.6).
Получающаяся в результате разностная схема является однородной и записывается на шаблоне типа крест. Соответствующая система разностных уравнений решается методом сопряженных градиентов Ortomin.
Для определения влияния неоднородности распределения температуры на характеристики получаемых эпитаксиальных слоев была проведена серия расчетов, в которых скорость охлаждения принимала значения a = 0.25 ,0.5 , толщина стенок составляла Ьь = 3mm. Параметры расплава приведены в таблице А.1.
Поскольку теплопроводность стенок достаточно велика, а толщина -мала по сравнению с горизонтальными размерами области, очевидно, что в процессе охлаждения, распределение температуры в подложке почти однородно, с небольшими возмущениями, локализованными около боковых стенок ростовой камеры. Величина этих отклонений слабо зависит от толщины жидкой фазы и линейных размеров области занятой расплавом и составляет при скорости охлаждения а = 0.25 5Т = 0.015, при а = 0.5 6Т — 0.03. Результаты расчетов показывают, что неоднородность температуры порождает пристеночный вихрь, имеющий форму вала, вытянутого вдоль боковой стенки. Затем движение распространятся внутрь области и приобретает форму пересекающихся валов, параллельных боковым стенкам (рис. В.21). Далее течение продолжает свое развитие, внутри области происходит смена планформы, а возле стенок структура течения остается достаточно устойчивой. Известно, что устойчивая структура течения вызывает сильную неоднородность формы поверхности выросшего слоя.
Из таблицы видно, что неоднородность, полученная с учетом конечной теплопроводности стенок и расплава, приблизительно в 1.5 раза выше, чем без ее учета. Для выяснения причин такого несоответствия были проведены расчеты со значениями чисел Gn, уменьшенными в полтора раза по сравнению с используемыми в п. 3.7. Оказалось, что влияние уменьшения Gn заметно только при небольшой толщине области.
1. Изучено влияние конвективного движения на процесс растворения подложки в методе жидкофазовой эпитаксии. Показано, что начальные возмущения в форме двумерных валиковых структур обладают большей скоростью роста, чем квадратные ячейки. При малой надкритичности (Ra 2.7Racr) это приводит к большей неоднородности поверхности травления. Поэтому они являются наиболее неблагоприятными с технологической точки зрения. При средней надкритичности (Ra « 4.3i?acr) неоднородность поверхности фронта травления практически не зависит от формы вносимого начального возмущения. Это позволяет при умеренном значении числа Ra 5Racr использовать двумерные модели для анализа задачи растворения подложки.
2. Наблюдаемые в расчетах существенно трехмерные структуры течения образованы ячейками g-типа, что полностью согласуется с данными теоретического анализа [71].
3. Показано, что конечная теплопроводность стенок формирует постоянно действующее возмущение, вызывающее существенный рост неоднородности подложки. Поэтому, при моделировании процесса эпитаксии с целью изучения однородности выросшего слоя, расчеты целесообразно проводить с учетом конечной теплопроводности стенок ростовой камеры.
