Моделирование роста кристаллов в условиях микрогравитации

-2-
СОДЕРЖАНИЕ
стр.
ВВЕДЕНИЕ 5
Г лава 1. Численный метод решения уравнений Навье-Стокса 21
§ 1. Формулировка проблемы 21
§ 2. Постановка задачи 26
§ 3. Преобразование операторов переноса 33
§ 4. Разнесенная сетка и способ расщепления уравнений 38
§ 5. Примеры пространственных аппроксимаций 46
§ 6. Окончательные разностные формулы, по которым проводился расчет задач 53
§ 7. Аппроксимация окончательными разностными формулами исходных уравнений 77
§ 8. Об устойчивости разностной схемы 89
§ 9. Обобщение численного алгоритма на трехмерный случай 95
Выводы к Главе 1 100
ГЛАВА 2. РЕШЕНИЕ ТЕСТОВЫХ ЗАДАЧ 102
§ 1. Стационарное течение несжимаемой жидкости 102
§ 2 Нестационарное течение сжимаемого газа 122
§ 3. Расчет течений сжимаемого газа. Сравнение численных схем 131
Выводы к Главе 2 138
-3-
ГЛАВА 3. МОДЕЛИРОВАНИЕ РОСТА ПОЛУПРОВОДНИКОВЫХ КРИСТАЛЛОВ. ЗАДАЧА О БЕСТИГЕЛЬНОЙ ЗОННОЙ ПЛАВКЕ КРЕМНИЯ 140
§ 1. Постановка задачи 140
§ 2. Продольные неоднородности и модели пограничного слоя 145
§ 3. Моделирование процесса бестигельной зонной плавки кремния на основе уравнений Навье-Стокса. Наземная отработка процесса 154
§ 4. Расчеты проведения процесса в условиях микрогравитации 169
Выводы к Г лаве 3 178
ГЛАВА 4. МОДЕЛИРОВАНИЕ РОСТА ПОЛУПРОВОДНИКОВЫХ КРИСТАЛЛОВ. ЗАДАЧА О ВЫРАЩИВАНИИ МОНОКРИСТАЛЛОВ СОЕДИНЕНИЙ А2В6 ФИЗИЧЕСКИМ ПАРОВЫМ ТРАНСПОРТОМ 180
§ 1. Постановка задачи 180
§ 2. Математическая модель процесса 185
§ 3. Моделирование земных условий роста. Рекомендации для проведения технологических процессов в условиях пониженной гравитации 192
§ 4. Граничные условия на поверхности пар- твердое тело 207
Выводы к главе 4 211
ГЛАВА 5. ЧИСЛЕННЫЙ МЕТОД ДЛЯ РАСЧЕТА ЗАДАЧИ РОСТА КРИСТАЛЛОВ В ДВУХФАЗНОЙ ОБЛАСТИ НА ПОДВИЖНЫХ СЕТКАХ 213
§ 1. Криволинейные координаты 213
§ 2. Решение задачи Стефана
225
-4-
§ 3. Задача о расчете температурного профиля электропечи “Кратер-ВМ”
Выводы к Главе 5
ГЛАВА 6. МОДЕЛИРОВАНИЕ РОСТА ПОЛУПРОВОДНИКОВЫХ КРИСТАЛЛОВ. ЗАДАЧА О НАПРАВЛЕННОЙ КРИСТАЛЛИЗАЦИИ АРСЕНИДА ГАЛЛИЯ В УСЛОВИЯХ МИКР01РАВИТАЦИИ
§ 1. Решение тепловой задачи. Поддержание постоянной скорости роста кристалла и управление формой фронта кристаллизации
§ 2. Экспериментальные результаты роста кристаллов методами направленной кристаллизации в условиях микрогравитации
§ 3. Решение задачи об осевой и радиальной микросегрегации при направленной кристаллизации нелегированного СаАэ в условиях пониженной гравитациии
Выводы к Главе 6
Заключение
230
242
244
244
265
278
295
298
Литература
305
-5-
ВВЕДЕНИЕ
Моделированию проблем тепломассообмена при росте кристаллов в космических условиях, т.е. решению задач гидромеханики, посвящено очень большое число работ. Имеется ряд монографий, как отечественных, так и зарубежных. Тем не менее, поток публикаций по этой теме не уменьшается, поскольку как сложность задач, так и требования к точности расчета непрерывно возрастают. Введение в эксплуатацию Международной космической станции значительно увеличило интерес к этим вопросам целого ряда экономически развитых стран.
Создание орбитальных станций и успешное осуществление длительных пилотируемых полетов привело в свое время к формированию нового научно-технического направления - космического материаловедения, которое поставило перед собой цель использовать состояние невесомости для получения веществ и материалов. В области материаловедения существуют направления, в которых использование космических условий наиболее перспективно. В первую очередь, это получение монокристаллов для интегральных схем, твердотельных лазеров и инфракрасной техники; материалов оптоэлектронной техники, сверхпроводящих материалов.
При постановке первых экспериментов космического материаловедения считалось, что на орбитальной станции можно легко вырастить однородные и совершенные кристаллы полупроводниковых материалов. Это уверенность основывалась на том предположении, что перенос массы и тепла в расплавах полупроводников будет происходить путем молекулярных процессов диффузии и теплопроводности, а тепломассоперенос путем конвекции будет полностью отсутс твовать или его вкладом можно будет пренебречь.
С начала проведения космических экспериментов было обнаружено, что при выращивании кристаллов, в первую очередь полупроводниковых
-6-
материалов, не достигается однородности электрофизических и структурных свойств. Как оказалось, на орбитальных станциях и других космических аппаратах режим невесомости не реализуется. Малые ускорения приводят к условиям, в которых конвекция в расплаве при росте кристаллов действительно становится малой, но, тем не менее, оказывает значительное влияние на свойства .выращиваемых кристаллов. Влияние микроускорений на рост кристаллов и их свойства является в настоящее время основным вопросом космического материаловедения. Новые знания, получаемые при проведении космических экспериментов в этой области, имеют фундаментальное значение и расширяют наши представления о сущности важнейших физических процессов.
Одним из первых открытий космического материаловедения явилось выяснение вклада термогравитационной конвекции (Марангони), т.е. конвекции негравитационного типа на рост и свойства кристаллов. Эта конвекция, обусловленная зависимостью градиента поверхостного натяжения свободной поверхности расплавов от температуры, теоретически была обнаружена уже давно, см. [1], но ее вклад считался пренебрежимо малым, а реальное влияние конвекции Марангони в земных условиях маскировалось действием термогравитационной и концентрационно-гравитационной конвекций. В космических условиях при незначительном влиянии сил гравитационной природы было обнаружено, что конвекция Марангони не уступает по влиянию гравитационной конвекции в наземных условиях.
Эксперименты, проводимые в космических условиях, требуют длительной подготовки и весьма дорогостоящи. Даже при активном освоении космического пространства, которое имело место в СССР и России за последние 25 лет, общее количество проведенных технологических экспериментов не так велико. Для полноценного анализа результатов, выяснения сущности явлений необходимо теоретическое исследование и обоснование проведения технологических экспериментов. Теоретическое
-7-
обоснование экспериментов играет в космическом материаловедении важнейшую роль и в большинстве случаев основывается на численном решении нелинейных уравнений в частных производных для расчета тепломассообмена в расплаве или газовой фазе при росте кристаллов.
Для первых лет развития космического материатоведения характерны относительно простые, полуаналитические модели процессов роста. В более позднее время анализ гидромеханических явлений и тепломассопереноса при росте кристаллов в условиях пониженной гравитации проводился, как правило, на основе решения уравнений Навье-Стокса в приближении Буссинеска, см. [2-5], а также [6-11]. Работы, представленные на последнем симпозиуме по гидромеханике гравитационно-чувствительных систем (Москва, ИПМех РАН, 2000г.), развивают и детализируют это направление, приводится решение некоторых трехмерных задач, см. [12-19]. В рамках решения этих уравнений проводилось численное исследование конвекции Марангони, [20-28].
Вместе с тем, был сформулирован ряд задач конвекции жидкости и газа, описание которых уравнениями Навье-Стокса в приближении Буссинеска непригодно. Здесь в первую очередь следует отметить работы М.Х. Стрельца,
А.И. Жмакина и Ю.Н. Макарова, см. [97-100], а также работы В.В. Пухначева, в последней из которых, см. [29], сформулирована трехмерная модель течения слабосжимаемой жидкости. В связи с этим возникла необходимость разработки численных схем для системы полных уравнений Навье- Стокса. Эта задача решается в рамках диссертационной работы (1985-1988г.г.). Разработана двумерная осесимметричная численная схема дозвуковых течений сжимаемого газа на основе полных уравнений.
Сложность технологических экспериментов, проводимых на орбитальной станции “Мир”, потребовали тщательного теоретического обоснования проведения наземной отработки и космических экспериментов. Некоторые из этих задач решались в диссертации с помощью уравнений Навье-Стокса для несжимаемой жидкости. Численному решению технологических задач
-8-
предшествовало подробное исследование точности метода на известных тестовых задачах. Основное внимание уделялось точности расчета течения внутри и вблизи пограничного слоя на неравномерных сетках, поскольку в задачах, связанных с ростом кристаллов, поведение расплава именно в этой области у фронта кристаллизации представляет основной интерес.
Еще одним фундаментальным результатом, наряду с обнаружением реального вклада конвекции Марангони и характеризующим особенности роста кристаллов в условиях микрогравитации, является обнаружение эффектов максимума температурной и концентрационной неоднородности, открытых
В.И. Полежаевым и сотрудниками, см.[30,31]. Эти эффекты получили название температурного (концентрационного) расслоения в замкнутых объемах. Их содержание заключается в наличии максимумов тепловой и концентрационной неоднородности при определенном значении гидродинамического числа подобия. Физически это означает, что механизмы конвекции и теплопроводности (диффузии) в технологическом процессе становятся одного порядка, и при этом достигается температурная (концентрационная) неоднородность большая, чем при преобладании любого из перечисленных механизмов. Как показала практика, максимум концентрационной неоднородности в большинстве космических экспериментов по росту полупроводниковых кристаллов .методами направленной кристаллизации, лежит как раз в диапазоне тех микроускорений (гидродинамических чисел подобия), при которых проводились эти эксперименты. Это обусловило значительно более высокое наличие радиальной макросегрегации в образцах указанных космических экспериментов, чем в их земных аналогах.
Теоретическое понимание этого фундаментального принципа позволило сформулировать к настоящему времени физическую картину роста кристаллов методами направленной кристаллизации в космических условиях и выяснить эффект влияния малых сил, поперечных направлению роста кристалла, B.C. Земсков и сотрудники, см. [32-35]. По мнению этих исследователей,
-9-
аналогичные эффекты могут иметь значение и для роста кристаллов в земных условиях при наличии малых сил, поперечных направлению роста кристалла.
Для адекватного теоретического обоснования проблемы о влиянии остаточных микроускорений, поперечных направлению роста кристаллов методами направленной кристаллизации, необходимо решать двухфазную задачу Стефана совместно с уравнениями, описывающими конвекцию в расплаве. Одной из первых работ по решению задачи Стефана является публикация А.А. Самарского и Б.Д. Моисеенко, см. [36]. В работе развит подход сквозного расчета уравнений с применением “размазывания” коэффициентов уравнений и, соответственно линии фронта кристаллизации, вдоль нескольких узлов сетки. Этот подход весьма удачен для решения задачи Стефана в области и сравнительно прост в реализации. Такой подход применяется и в настоящее время, см., например, работы А.С. Овчаровой [37,38]. Существенным недостатком этого подхода является, ввиду отсутствия выделенного фронта кристаллизации, невозможность рассчитать течения в узком пограничном слое у фронта кристаллизации и, таким образом, решить задачу о макросегрегации или микросегрсгации примеси.
Представляет интерес подход, развитый С. Lan, см. [39-41], при решении двумерной, а потом трехмерной задачи Стефана для бестигельной зонной плавки. Развивается более детальное описание, чем простое “размазывание” коэффициентов и линии фронта кристаллизации. Используется методика “tracking of crystallization front”, при котором различаются области расплава и кристалла, но линия фронта кристаллизации также явно не выделяется. При моделировании трехмерной задачи вертикального роста кристаллов решалась нестационарная задача движения расплава, переноса тепла и движения поверхности раздела. Обсуждались вопросы, связанные с бифуркацией двумерных решений при переходе к трехмерным в случае возрастания интенсивности конвекции.
-10-
Из работ, в которых явно выделяется фронт кристаллизации, следует отметить решение трехмерной задачи Стефана при замерзании воды, см. [42]. Задача решалась в переменных вектор вихрь - векторный потенциал, используя разработанный V. Davis et al. код FREEZE3D.
Среди отечественных работ, используемых при моделировании технологических процессов, известен пакет двумерных программ КАРМА 1 (И.В. Фрязинов и соавторы), см. [43], где приводится нестационарное решение полной системы уравнений. Комплекс базируется на неявной схеме без расщепления в переменных функция тока - вихрь. Для решения возникающей алгебраической системы уравнений используется прямой метод. В диссертации применяется другой подход, разработанный примерно в одно время с [43]. Он основан на применении схемы расщепления в естественных переменных, которая обобщена на случай подвижных непрямоугольных сеток. Специфика численного метода, построенного в диссертации, позволила разработать принципиально новый способ решения задачи Стефана, при котором одновременно вычисляются изменения теплового поля за шаг по времени и новая скорость роста кристалла.
Предлагаемый в диссертации метод решения задачи Стефана позволяет заранее задавать способ изменения внешнего теплового поля (температурных граничных условий) для поддержания в процессе постоянной наперед заданной скорости роста кристалла. Этот метод решения двухфазной задачи дает возможность рассчитать как радиальное и осевое отклонение примеси (макросегрегация), так и численно исследовать колебания скорости роста кристалла во времени (полосы роста, или микросегрегация). Такие возможности построенного численного метода решения задачи Стефана являются необходимыми при анализе роста кристаллов методами направленной кристаллизации в космических условиях.
В целом, основной проблемой, стоящей перед космическим материаловедением в течение продолжительного ряда лет, является
-11-
исследование влияния остаточных микроускорений и связанного с ними гидродинамического фактора на рост и свойства кристаллов, выращенных в космических условиях. Диссертация посвящена решению этой проблемы. Разработаны современные численные методы расчета влияния конвекции на свойства выращенных кристаллов, проведено теоретическое исследование и обоснование технологических экспериментов на борту орбитальной станции “Мир” и сопоставление их с другими космическими экспериментами по росту кристаллов; выработаны рекомендации для проведения наземных процессов и экспериментов на Международной космической станции.
Целью диссертационной работы являлось развитие научного направления космического материаловедения, связанное с теоретическим обоснованием и численным исследованием ряда актуальных задач по росту кристаллов в космических условиях, в том числе:
-исследование влияния гидродинамических факторов на рост и свойства выращенных кристаллов при наземном росте и в условиях космических экспериментов путем численного моделирования технологических процессов; исследование гравитационной чувствительности полупроводниковых и ряда других кристаллов;
-разработка современных численных методов решения задач конвекции на основе полных уравнений Навье-Стокса для сжимаемого газа и в приближении Буссинеска для несжимаемой жидкости при росте кристаллов; разработка методов совместного решения задачи Стефана для двухфазных сред и задачи конвекции в расплаве с подвижным фронтом кристаллизации;
-теоретическое обоснование технологических экспериментов на орбитальной станции “Мир”, сопоставление их с другими космическими экспериментами по росту кристаллов; выработка рекомендаций для проведения экспериментов на Международной космической станции.
-12-
Диссертация состоит из введения, шести глав, заключения и списка использованных источников. Объем диссертации 326 страниц, включая 52 рисунка, 15 таблиц и библиографию из 187 источников. В тексте принята следующая нумерация: в формулах первая цифра означает номер главы, вторая-номер параграфа, третья - является порядковым номером формулы в параграфе. Для рисунков и таблиц применяется двойная нумерация: первая цифра- номер главы, вторая- номер рисунка или таблицы в главе, соответственно.
Первая глава диссертации посвящена разработке численного метода решения задач космического материаловедения, во второй главе проводится тестирование разработанной численной схемы. С помощью этого метода в главе 3 и 4 решены задачи о бестигельной зонной плавке кремния и физическом паровом транспорте соединений А^В6 в парах собственного компонента. Для решения задачи о нормальной направленной кристаллизации арсенида галлия в главе 6 необходимо было предварительно разработать и описать в главе 5 метод решения задачи Стефана и рассчитать тепловые параметры электропечи. В заключении приведены общие выводы диссертационной работы.
Несмотря на достаточно развитые методы решения уравнений Навье-Стокса, большинство задач расчета конвективных течений решается в приближении несжимаемой жидкости. В первой главе подробно описан численный метод решения полных уравнений Навье-Стокса для сжимаемого газа, сопоставимый по эффективности с численными схемами для уравнений несжимаемой жидкости в приближении Буссинеска. Этот численный метод можно успешно применять как для расчета газодинамических течений в космических условиях, так и для решения обширного класса наземных задач дозвукового течения газа. Приведена численная схема, в которой выполняются закон сохранения массы, балансы кинетической и внутренней энергии газа. Производится расщепление уравнений по физическим процессам, а там, где это оправдано, и по пространственным переменным. На центральном этапе
-13-
расщепления предложенной схемы организован эффективный итерационный процесс, на каждом шаге которого необходимо решать уравнение типа Пуассона. Численный метод для уравнений сжимаемого газа построен с учетом специфики дозвуковых течений и сопоставим по эффективности численным схемам для уравнений несжимаемой жидкости.
Разработанный численный алгоритм для решения полных уравнений Навье-Стокса построен таким образом, что вычислительная программа может быть легко переключена на расчет уравнений несжимаемой жидкости. Расчет сильносжимаемых течений предложенным методом занимает всего в 2-2,5 раза времени больше, чем расчет в рамках приближения Буссинеска. Недостатком ряда численных схем в приближении Буссинеска является сложность постановки граничных условий на твердой поверхности, каким является, например, фронт кристаллизации. В построенном численном методе эта трудность отсутствует, что является особенно важным не только при моделировании конвекции в условиях пониженной гравитации, но и при проведении отработки технологических процессов в земных условиях, которые характеризуются быстрым течением расплава в тонких пограничных слоях.
Вторая глава посвящена тщательному тестированию построенного численного метода как на задачах с существенно сжимаемыми течениями, так и на задачах для несжимаемой жидкости. При решении ряда тестовых задач, приведенных в этой главе, показана высокая точность расчета задач конвекции предложенным численным методом и разработаны приемы решения задач в областях с тонкими пограничными слоями на неравномерных сетках.
В третьей главе решена задача о бестигельной зонной плавке кремния в земных условиях и в условиях пониженной гравитации. Задача решалась наряду с наземной отработкой этого технологического процесса на земном аналоге аппаратуры “Оптизон” со световым нагревом зоны расплава и космическими экспериментами, выполненными на этой аппаратуре, расположенной на борту орбитальной станции “Мир”.
-14-
При решении задачи в земных условиях исследовались сложные механизмы распределения примеси в зоне расплава и в тонком пограничном слое у фронта кристаллизации, которые определяются результатом взаимодействия конвективной ячейки и диффузионных процессов. В рассмотренном режиме тепломассопереноса важную роль в распределении примеси играют апериодические колебания конвективной ячейки около своего положения равновесия. В пограничном слое у фронта кристаллизации подробно описана радиальная неоднородность примеси, которая приводит на практике к неоднородности электрофизических и структурных свойств выращенных кристаллов.
При решении задачи в условиях пониженной гравитации обнаружено, что режимы микроускорсний, имевшие обычно место на орбитальной станции “Мир”, соответствуют максимуму концентрационного расслоения. В результате у фронта кристаллизации создается значительно большая радиальная неоднородность примеси, чем в земных условиях. Используя теорию подобия, сделаны выводы о характере роста в условиях микрогравитации кристаллов большого промышленного диаметра.
При решении задачи о бестигельной зонной плавке кремния условия получения диффузионного режима роста были сформулированы в терминах конкретных значений гидродинамических параметров подобия (Рейнольдса или, соответственно, Рэлея). В настоящее время эти значения числа Рэлея являются общепризнанными критериями успешного роста полупроводниковых кристаллов в условиях микрогравитации.
Четвертая глава посвящена решению задачи о росте кристаллов соединений А2В6 физическим паровым транспортом в парах собственного компонента на монокристаллическую затравку. При решении этой задачи сформулирована физическая картина явления, а расчеты по численной модели позволили оптимизировать параметры технологического процесса. Обнаружено, что в отсутствии инертного газа в ампуле основным содержанием
-15-
процесса является непосредственный перенос пара с источника на растущий кристалл под воздействием возникающей на них разности давлений собственных паров. Этот перенос определяется граничными условиями на поверхности пар - твердое тело, следующими из кинетики межфазных процессов, а также общим давлением в системе, которое поддерживается величиной температуры в низкотемпературной части ампулы.
При решении задачи получено, что в случае большого общего давления в ампуле реализуется диффузионный режим роста, но скорость роста при этом невелика и поэтому не удается вырастить объемные кристаллы значительных размеров. При уменьшении общего давления в ампуле благодаря непосредственному массопереносу с источника на кристалл, обусловленного заметной разностью давлений, скорость роста резко возрастает. Но при этом значительно ухудшается качество кристаллов. Выяснено, что практическая задача заключается в выборе промежуточного режима, когда при умеренной скорости роста удается выращивать высококачественные объемные монокристаллы с высоким удельным сопротивлением.
В задаче оптимизированы режимы роста кристаллов селенида кадмия на аппаратуре “Кратер”. Вычислены значения температур, при которых, как показано в технологических экспериментах, выращиваются близкие к стехиометричным высокоомные объемные монокристаллы с достаточно высокой скоростью роста. Кроме того, решение задачи позволило выяснить, что при росте кристаллов этим методом непосредственный перенос вещества с источника на кристалл не имеет гравитационной природы. Поэтому были сделаны выводы о нецелесообразности проведения подобных технологических процессов в условиях пониженной гравитации.
Значительную часть космических экспериментов, как по общему количеству, так и по их важности для космического материаловедения, составили процессы роста кристаллов полупроводниковых соединений методами направленной кристаллизации. Для выяснения причин значительной
-16-
примесной неоднородности в этих экспериментах и подтверждения физической картины макросегрегационных процессов, разработанной к настоящему времени, необходимо решать двухфазную задачу Стефана. В пятой главе описан новый метод решения задачи Стефана на подвижных сетках, который включает в себя решение системы уравнений Навье-Стокса в расплаве, уравнения теплопроводности в кристалле, условия Стефана и оттеснения примеси на заранее неизвестном неплоском фронте кристаллизации. Специфика численного метода решения уравнений Навье-Стокса в расплаве позволила разработать такой метод решения задачи Стефана, при котором одновременно вычисляются изменения теплового поля за шаг по времени и новая скорость роста кристалла.
Проводится обобщение численной схемы расчета конвекции в расплаве, разработанной в первой главе, на случай непрямоугольных сеток. Введены криволинейные координаты и контравариантные компоненты скоростей. При построении сетки существенно используется особенность задач роста полупроводниковых и некоторых других кристаллов, у которых искривление фронта кристаллизации не может быть слишком велико. В результате разностная схема использует сильно неравномерную сетку при сравнительно простом способе ее построения и небольшом общем числе узлов.
В третьем параграфе выполнена оптимизация тепловых режимов работы электропечи, предназначенной для роста кристаллов в космических условиях. Рассчитан теплообмен и построена математическая модель для управления тепловым полем электропечи “Кратер-ВМ”. Вычислены коэффициенты теплоотдачи каждого нагревателя, рассчитана плотность потока тепла через пространство между нагревателями. Это позволило реализовать необходимые температурные профили в высокотемпературной части электропечи, рассчитать величину необходимого расстояния между высокотемпературной и низкотемпературной частями электропечи.
-17-
Шестая глава посвящена решению задачи Стефана при росте кристаллов арсеиида галлия на аппаратуре “Кратер-ВМ”. В первом параграфе рассчитаны тепловые аспекты задачи при ее наземной отработке, найдены способы поддержания постоянной скорости роста кристалла и управления формой фронта кристаллизации, что регулирует длину монокристаллической части слитка. Проводится сравнение расчетных данных с непосредственными экспериментами по росту кристаллов, приводятся характеристики кристаллов арсенида галлия, полученные в ходе проведения технологических процессов.
Во втором параграфе главы перечислены некоторые космические эксперименты по росту кристаллов методами направленной кристаллизации и приведена физическая модель макросегрегации, разработанная Российскими исследователями, которая объясняет причины значительной радиальной неоднородности примеси в этих космических экспериментах. Приводится основной принцип, которому подчиняется примесная неоднородность в расплавах - принцип максимума концентрационной неоднородности в замкнутых объемах. Подтверждено, что большинство проведенных космических экспериментов лежит в диапазоне чисел подобия, отвечающих неблагоприятному режиму принципа концентрационной неоднородности. Приведена оценка чисел подобия, совпадающая по порядку величины с расчетными результатами в третьей главе диссертации, при которой массоперенос в расплаве становится чисто диффузионным.
В третьем параграфе главы сформулирована и решена задача Стефана для исследования микросегрегации (полос роста) при направленной кристаллизации арсенида галлия в условиях орбитальной станции “Мир”. Эта задача косвенно подтвердила основные положения физической модели процесса, сформулированной другими авторами для условий макросегрегации. Кроме того, были исследованы процессы, связанные с колебаниями скорости роста, а также осевой и радиальной сегрегацией, которые являются следствием этих колебаний.
-18
На защиту выносятся:
1. Численный метод расчета конвекции в сжимаемом газе на основе решения полной системы уравнений Навье-Стокса;
численный метод решения уравнений несжимаемой жидкости в приближении Буссинеска, позволяющий получать высокую точность при расчете конвекции в тонких пограничных слоях.
2. Закономерности распределения примеси при бестигельной зонной плавке кремния в зоне расплава и особенности радиальной сегрегации в тонком пограничном слое у фронта кристаллизации;
получены зависимости максимума концентрационного расслоения в условиях пониженной гравитации от значения гидродинамического числа подобия и определены условия диффузионного режима роста для кристаллов различного диаметра.
2 6
3. Механизмы роста кристаллов соединений А В при физическом паровом транспорте в парах собственного компонента на монокристалл и ческу ю затравку;
оптимизирован рост кристаллов селенида кадмия. Вычислены технологические параметры, при которых выращиваются близкие к стехиометричным высокоомные объемные монокристаллы с достаточно высокой скоростью роста.
4. Метод совместного решения задачи Стефана для двухфазных сред и задачи конвекции в расплаве с подвижным фронтом кристаллизации;
проведена оптимизация температурного режима и управления тепловым полем аппаратуры, предназначенной для роста кристаллов в условиях микрогравитации.
5. Условия поддержания постоянной скорости роста кристалла и способ управления формой фронта кристаллизации при росте кристаллов методом
-19-
направленной кристаллизации в земных условиях и в условиях пониженной гравитации;
механизм и условия образования продольных и поперечных полос роста при выращивании кристаллов арсенида галлия в условиях остаточных вибрационных микроускорений.
Практическая значимость полученных результатов:
1. Разработан численный метод и соответствующий ему комплекс программ для моделирования дозвуковых течений сжимаемого газа на основе полной системы уравнений Навье-Стокса, сравнимый по эффективности с численными схемами для решения задач конвекции несжимаемой жидкости.
2. Определены механизмы переноса примеси при бестигельной зонной плавке кремния, позволившие сформулировать физическую картину процесса и сформулировать микрогравитационные условия, при которых будет достигаться рост однородных и совершенных кристаллов.
3. Выяснена физическая картина роста кристаллов соединений А2В6 физическим паровым транспортом в парах собственного компонента, что дало возможность провести оптимизацию земных процессов роста и показать неперспективность проведения подобных ростовых процессов в условиях микрогравитации.
4. Разработан совместный метод решения задачи Стефана в двухфазной области и задачи конвекции в расплаве, позволяющий эффективно численно решать и оптимизировать задачи роста полупроводниковых и ряда других кристаллов методом направленной кристаллизации; управлять технологическим процессом роста в условиях микрогравитации.
5. В результате решения задачи Стефана подтверждена выработанная к настоящему времени физическая модель макросегрегационных явлений при росте кристаллов методами направленной кристаллизации в космических условиях, которая объясняет причины значительной радиальной
-20-
неоднородности примеси в большинстве проведенных этими методами космических экспериментов. Показано, что характер образования полос роста в поле низкочастотных вибрационных микроускорений аналогичен указанным макросегрегационным явлениям.
Апробация работы. Результаты диссертации докладывались и обсуждались на: IV Всесоюзной конференции “'Гермодинамика и материаловедение полупроводников” (Москва, 1989); III Всесоюзной конференции “Моделирование роста кристаллов” (Рига, 1990); I Международном симпозиуме по гидромеханике и тепломассообмену в микрогравитации (Пермь- Москва, 1991); I Аэрокосмической конференции “Перспективы производства в космосе” (Москва, 1992); I Международном аэрокосмическом конгрессе (Москва, 1994); I Международной авиакосмической конференции “Человек- Земля- Космос” (Москва, 1995); конференции Отделения Микроэлектроники и Информатики Международной Академии Информатизации (Москва, 1997); Объединенном X Европейском и VI Российском симпозиуме “Физические науки в микрогравитации” (Ст.-Петсрбург, 1997); Научно-исследовательском семинаре “Механика невесомости и гравитационно- чувствительные системы” (Москва, 1998); VII Российском симпозиуме “Механика невесомости. Итоги и перспективы фундаментальных исследований гравитационно- чувствительных систем” (Москва, 2000); IX Национальной конференции по росту кристаллов НКРК-2000 (Москва, 2000); III Международной научно-технической конференции “Электроника и информатика- XXI век” (Москва, 2000).
Основные результаты работы опубликованы в 21 статье в реферируемых журналах, 1 научно-техническом отчете и доложены в 13 выступлениях на конференциях и симпозиумах [44-78].
-21-
Глава 1. Численный метод решения уравнений Навье-Стокса § 1. Формулировка проблемы
Явление естественной конвекции при росте кристаллов из расплава или газовой фазы описывается общей системой уравнений газовой динамики. Эта система включает в себя уравнения Навье-Стокса. уравнение баланса энергии, уравнение диффузии примеси, уравнение непрерывности и уравнение состояния [79,80]:
д v 4
р (— + (V • V)v) = - grad р + pg + p-(- graddivv - rotrot v), (, д i)
3 T
pCv(— + V’VT) = k-divgradT-P’divv+ fiD^ ^ { 2)
dt
o(-
dt
3C
P (— + v • VC) = D • div (p gradC), (1 л
f + d;v(Pv) = 0, (U4)
p = p(P,T). (1.1.5)
Здесь V - вектор скорости, р - давление, Т - температура, р -плотность, § - ускорения свободного падения, р- коэффициент
динамической вязкости, Су- удельная теплоемкость при постоянном объеме, к - коэффициент теплопроводности, О- коэффициент диффузии, С- концентрация примеси, отвечает за работу сил вязкости.
Эти уравнения обладают большой степенью общности и описывают широкий класс течений газа (жидкости). Уравнения (1.1.1)-(1.1.5) достаточно сложны, и до настоящего времени для исследования конвекции обычно применяют упрощенные уравнения, а именно систему уравнений Навье-Стокса в приближении Буссинеска.
Основная особенность этого приближения заключается в учете сжимаемости среды только в массовой силе уравнения (1.1.1), а во всех
-22-
других уравнениях плотность среды полагается постоянной. Несмотря на кажущуюся непоследовательность, это приближение очень хорошо себя зарекомендовало, и большая часть вычислительной информации о естественной конвекции за последние десятилетия получена с помощью решения уравнений в приближении Буссинеска [81-94,43]. Остановимся на ограничениях приближения Буссинеска более подробно.
Во-первых, для заметно сжимаемых сред (или достаточно больших градиентов температуры в области) приближение, вообще говоря, неприменимо. Часто расчеты в приближении Буссинеска проводят и для этих случаев, но достоверность получаемых результатов остается спорной.
Во-вторых, при решении уравнений в приближении Буссинеска существуют определенные вычислительные трудности. Здесь возможны два основных подхода, связанные с выбором расчетных переменных. При решении уравнений в естественных переменных “скорость-давление” основным преимуществом является простота постановки граничных условий для поля скорости как в случае твердых стенок, так и на свободных поверхностях. Существенной трудностью при использовании естественных переменных является необходимость решать уравнение Пуассона с условиями Неймана на всех границах, что является сложной вычислительной задачей (см., например, [84,93]). Для другого подхода, использующего переменные “вихрь- функция тока”, указанную вычислительную трудность удается избежать. Для “функции тока” также необходимо решать уравнение Пуассона, но уже с граничными условиями Дирихле, что является значительно более легкой вычислительной задачей. В то же время такой подход имеет существенный недостаток, связанный с появлением дополнительного граничного условия для функции “вихрь”, которого нет в исходной постановке задачи, см. [90,93,84]. Как показали приведенные в главе 2 результаты сравнения решения тестовой задачи с расчетами других авторов, это приводит к заметной погрешности
-23-
вычисления решения вблизи твердой стенки. Наконец, возможен полностью неявный метод решения задачи с одновременным вычислением переменных “вихрь” и “функция тока”, см. [43]; но это является громоздкой вычислительной задачей и применяется редко.
В целом, несмотря на значительные успехи, существует настоятельная необходимость в разработке численных методов как для системы уравнений газовой динамики, т.е. полной системы уравнений Навье-Стокса (1.1.1)-(1.1 -5), так и для решения уравнений в приближении Буссинеска. Рассмотрим сначала различные методы решения полной системы уравнений.
Явная схема [95] для решения полной системы уравнений Навьс-Стокса проста и экономична. Но условие ее устойчивости ограничивает шаг по времени величиной т < Н/с, где 1г- шаг по пространственной координате, с- скорость звука. Характерное время самого конвективного процесса значительно больше и равно т = Н ! V, где V - гидродинамическая скорость газа. Таким образом, явная схема [95] для расчета конвективных процессов в сжимаемом газе неэффективна.
Естественно потребовать, чтобы разностная схема для полных уравнений Навье-Стокса была бы работоспособной и в предельном случае практически несжимаемого течения, не слишком уступая по эффективности разностным схемам для уравнений в приближении Буссинеска. Это означает: I) должен быть сравним объем вычислительной работы при расчете одного шага по времени; 2) расчеты с одним и тем же шагом по времени должны обеспечивать одинаковую точность.
В случае разностных схем для уравнений в приближении Буссинеска каждый шаг по времени связан, как правило, с решением разностного аналога уравнения Пуассона. Для уравнений сжимаемого газа существует схема полного расщепления по пространственным переменным [96]. Эта схема на одном шаге по времени даже более эффективна, чем схемы для
-24-
уравнений несжимаемой жидкости. Но второе условие оказывается невыполненным. Это связано с необходимостью проведения расчетов с шагом по времени т « hi v. Такой шаг отвечает большим значениям числа Куранта (К = г - с/ И\ что приводит к увеличению той части погрешности аппроксимации, которая вводится в схему процедурой расщепления. Это было наглядно продемонстрировано в [96], когда расчет с шагом по времени, г«4г0, где г0- шаг по времени, который следует из условия Куранта, приводил к искажению решения. Авторами [96] была предложена трехслойная схема переменных направлений, для которой эффект понижения точности отсутствовал вплоть до г»8г0. Тем не менее при т « 16г0 точность расчета опять заметно падала.
Подробнее, причина потери точности заключается в следующем. В случае дозвукового слабосжимаемого течения в уравнении непрерывности
—+ v • grad р) = -div v
р dt
левая часть много меньше, чем каждая из пространственных производных, входящих в div У. В результате при т » т0 погрешность аппроксимации в
основном определяется той своей частью, которая вносится в схему процедурой расщепления но направлениям. Поэтому схемы для дозвуковых течений следует строить аналогично схемам для несжимаемой жидкости: оба слагаемых, входящих в div У, следует аппроксимировать на одном временном слое. Хотя на каждом временном слое приходится решать двумерную систему алгебраических уравнений, допустимый шаг по времени практически перестает зависеть от условия Куранта, как это будет показано в главе 2.
Трудности в построении эффективной схемы для дозвуковых течений сжимаемого та в рамках полных уравнений Навье-Стокса привели к созданию приближенной теории [97-100], где разностные схемы
-25-
с нужными свойствами строить проще. Имеются несомненные успехи в применении этой теории к некоторым задачам. Тем не менее, для получения точного решения возникающих в этой области сложных и разнообразных задач конвекции необходима разработка разностных схем расчета полных уравнений сжимаемого газа.
В первой главе диссертации формулируется численный метод исследования конвекции в газе, который основан на решении системы полных уравнений (1.1.1)-(1.1.5), то есть плотность считается переменной во всех уравнениях системы. В качестве уравнения состояния р = р{ТуР) выбрано уравнение состояния идеального газа
п _ К т
7/’ (1.1.6)
где Я- универсальная газовая постоянная, р0- молекулярный вес.
Разработан эффективный метод решения системы уравнений (1.1.1)-(1.1.4),
(1.1.6). Как показало сравнение, расчет сильносжимаемых течений занимает всего в 2-2,5 раза времени больше, чем аналогичный расчет в рамках приближении Буссинеска.
Значительным достоинством построенного численного метода является тот факт, что он включает в себя возможность расчета несжимаемых задач. Для этого достаточно один из этапов расщепления задачи заменить на другой, более простой. В случае приближения Буссинеска разработанный алгоритм использует удачные стороны двух способов расчета конвекции. Используются естественные переменные, однако вместо уравнения Пуассона с условиями Неймана решается уравнение Пуассона с условиями Дирихле для разностного аналога функции тока.
-26-
Численному решению задач о росте полупроводниковых кристаллов предшествовало тщательное исследование точности метода на известном международном тесте (глава 2, §1). Основное внимание уделялось точности расчета течения внутри и вблизи пограничного слоя на неравномерных сетках, поскольку в прикладных задачах поведение расплава именно в этой области у фронта кристаллизации представляет основной интерес. Как показали расчеты, предложенный алгоритм решения задачи на сравнительно небольшом числе узлов дал значительно более точные результаты, чем решение задачи другими численными методами.
Запишем размерные уравнения естественной конвекции сжимаемого
газа:
§ 2. Постановка задачи
р С,, (— + V • V Г) = А • сЦу%гаАТ - р • сИх р О ь
Л# *
0/
(1.2.2)
р( + у-УС) = й-с/м(р grad С),
д!
(1.2.3)
0Г
(1.2.4)
Р = ~рТ.
(1.2.5)
-27-
Здесь V- вектор скорости с компонентами (и,\у) вдоль осей ег,ег соответственно (цилиндрические координаты), через обозначено
выражение
гл п//^мч2 /3 и\2 и2. .ди ди\2 2,1 дги дм\2
зГ' + ^ *& + аГ* "з 7"аГ+ з7 <>-2«
что соответствует работе сил вязкости. Переменные /' и 2 в формуле
(2.1.6)- цилиндрические координаты.
Численный метод исследования конвекции в газе, предлагаемый в диссертации, основан на решении системы полных уравнений (2.1.1)-
(2.1.5), то есть плотность считается переменной во всех уравнениях системы. Предлагаемый численный метод [44-46,65] основан на расщеплении уравнений как по физическим процессам, так и, где это оправдано, по пространственным переменным. Для получения второго порядка точности по времени применяется схема симметричного расщепления [101]. Среди физических процессов выделены задача о переносе вдоль траектории, задача о вязком обмене и так называемая задача о релаксации, в которой из пространственных слагаемых в уравнении движения участвуют градиент давления и массовая сила.
После расщепления по пространственным переменным задачи о переносе вдоль траектории и о вязком обмене реализуются с помощью одномерных прогонок.
На этапе о вязком обмене дополнительно выделены слагаемые, возникающие только в сжимаемом случае. Наличие смешанных производных не позволяет применить для их расчета алгоритм прогонки. Был организован итерационный процесс, являющийся модификацией метода переменных направлений. Для сходимости обычно достаточно не более трех итераций.
На этапе расщепления задачи о релаксации из уравнений движения участвуют градиент давления и массовая сила, из уравнения баланса
-28-
энергии - работа сил давления и - слагаемые, уравнения непрерывности
и состояния присутствуют полностью.
Все изменения плотности в уравнениях происходят на этапе задачи о релаксации. Таким образом, организация вычислительного алгоритма для получающейся на этом этапе системы алгебраических уравнений является узловым моментом построения численного метода. Организован оригинальный итерационный процесс, использующий процедуру исключения давления, характерную для уравнений несжимаемой жидкости.
Вычислительный алгоритм на этапе задачи о релаксации организован следующим образом. Искомыми здесь являются разностные функции V, р, Т и р. Воспользуемся известным приемом [102], когда уравнения разбиваются на две группы. В первую входит расщепленное уравнение баланса энергии, а во вторую- оставшиеся уравнения. Сначала при заданной функции Т решаются уравнения второй группы и находятся остальные неизвестные функции. Затем из уравнения баланса энергии находится Т и так далее. При проведении расчетов задавалось максимально допустимое относительное изменение температуры на одной
внешней итерации д -10“3, для достижения которой оказалось достаточным одной-двух внешних итераций. Эффективность численного метода в основном определяется итерационным процессом для решения второй группы уравнений.
Итерационный процесс (по плотности) для второй группы уравнений строится в два этапа. Сначала из расщепленных уравнений движения и уравнения непрерывности исключаются давление Р и компонента и скорости V, а для компоненты XV скорости V получаем разностное уравнение типа Пуассона. Плотность р берется с нижнего слоя. В случае равномерных сеток уравнение можно решить прямым методом.
-29-
Второй этап итерации по плотности состоит в следующем. Давление Р из уравнения состояния подставляется в расщепленные уравнения движения. Плотность р выражается из уравнения непрерывности. В результате получается система двух линейных уравнений относительно и и \у. Из этих уравнений получаются поправки к предыдущим значениям и и \у. По известным скоростям и и \у вычисляется из уравнения непрерывности плотность р. На этом итерация по плотности закончена, вычислены новые значения и, \у и р.
Во всех проведенных расчетах этот итерационный процесс сходился достаточно быстро, за 2-3 итерации, даже при сравнительно больших (0.8-1.0) числах Маха. Учитывая внешние итерации, связанные с изменением температуры, на одном шаге по времени приходилось не более шести раз решать уравнение Пуассона. Как показало сравнение, расчет сильносжимаемых течений занимает всего в 2-2,5 раза времени больше, чем аналогичный расчет в рамках приближения Буссинеска.
Большим достоинством построенного численного метода является тот факт, что он включает в себя возможность расчета несжимаемых задач (приближение Буссинеска). Для этого достаточно этап расщепления задачи о релаксации заменить на другой, более простой. В несжимаемом случае уравнение неразрывности позволяет ввести на этом этапе расщепления разностный аналог функции тока ц/. Исключая давление Р из расщепленных уравнений движения, для определения функции тока ц/ получается разностное уравнение типа Пуассона с граничными условиями Дирихле. Таким образом, в несжимаемом случае на слое по времени достаточно один раз обратить оператор Лапласа.
Возможность переключения управляющей программы на расчет полных уравнений или уравнений в приближении Буссинеска представляет большое удобство.
-30-
Во-первых, для несжимаемого варианта (программы уравнений в приближении Буссинеска) существует очень удобный стационарный тест (У.Оау1з) для количественной проверки точности расчетов по численной схеме как для наземных условий проведения технологических процессов, так и для условий микрогравитации.
Во-вторых, несжимаемый вариант позволяет проверить разностную схему на уже исследованных сложных двумерных нестационарных задачах.
В-третьих, при решении сжимаемого варианта (численной схемы для полных уравнений) часто бывает важно знать, как и насколько отличается сжимаемый вариант от несжимаемого. Параллельный расчет обоих вариантов дает ответ на этот вопрос.
Все эти возможности, перечисленные выше, были реализованы в диссергационной работе.
Для численного решения уравнений (1.2.1)-(1.2.5) введем безразмерные переменные (они обозначены штрихами) по следующим формулам: г = 2г', z = Zz,, где Z~ характерный размер области, v = Uv\t = Zt'IU. В качестве характерной скорости выбрано 11 = рgZ, где /3 = АТ/Т[- относительный перепад температуры в области. Для плотности, температуры, концентрации примеси и давления получим р — р0 р\ Т = АТ • V + 7], С = АС' + С0, р = Р р'. Здесь р0-характерная плотность газа, АТ и АС- перепад температуры и концентрации в области, 7^ и С0- размерные значения температуры и концентрации примеси, Р = /?01! .
При нормальных условиях (плотности р0 и температуре Т{) безразмерное давление в уравнении состояния равно