Содержание к диссертации
Введение
2. Постановка задачи 34
2.1. Математическая модель процессов гидродинамики и теплообмена на основе уравнений Навье-Стокса (приближение Буссинеска) 34
2.2. Математические модели переноса примеси 40
2.3. Критериальное уравнение и диапазон изменения параметров 43
3. Методы численного моделирования 47
3.1. Метод конечных разностей 48
3.2. Метод конечных элементов 60
3.3. Некоторые сведения о программной реализации методов 69
3.4. Тесты численных решений 71
4. Изотермические течения 84
4.1. Вращение кристалла 84
4.2. Противовращение кристалла и тигля 94
4.3. Особенности течения в двойном тигле 99
4.4. Влияние изотермического течения на распределение примеси 101
5. Неизотермические течения 106
5.1. Тепловая конвекция 106
5.2. Некоторые эффекты совместного действия тепловой конвекции и вращения кристалла 112
5.3. Особенности течения и теплообмена при совместном действии тепловой' конвекции и противовращения кристалла и тигля 122
5.4. Влияние неизотермичности течения на распределение примеси 134
5.5. О выборе параметров для выращивания кристаллов методом Чохральского 140
Заключение 143
Литература
- Математические модели переноса примеси
- Метод конечных элементов
- Противовращение кристалла и тигля
- Некоторые эффекты совместного действия тепловой конвекции и вращения кристалла
Введение к работе
Предмет исследования
Совершенные кристаллы находят все большее применение в технике, будучи основой для создания интегральных схем ЭВМ, лазерных установок, фотоприемников и средств связи. Одновременно повышаются требования к совершенству кристаллов и увеличивается объем их производства. Большинство кристаллов выращивается из расплава методами Бриджмена, бестигельной зонной плавки и методом Чохральского. Последний является самым распространенным промышленным методом получения кристаллов полупроводниковых, оптических, драгоценных и многих других кристаллов.
Совершенство кристаллов зависит от состава расплавленной массы и процессов тепло- и массообмена на стадии предшествующей кристаллизации (см. Лодиз, Паркер [i] ). Кристаллы различаются по физическим свойствам и температуре кристаллизации, что приводит к отличиям в установках и параметрах выращивания различных материалов одним и тем же методом, например, методом Чохральского. Одной из важных особенностей процессов роста для различных материалов является характер подвода тепла, его величина и ориентация по отношению к фронту кристаллизации и силе тяжести, что приводит к существенно различным структурам движения расплава, интенсивности тепло- и массообмена и, в свою очередь, оказывает различное влияние на макро- и микроструктуру кристаллов.
Вопросы термодинамики и фазового равновесия при росте кристаллов, а также влияние на процесс кристаллизации движения расплава рассмотрены в книге Ro$ene.m2i a\z\ . Обзоры возможных механизмов движения при выращивании кристаллов из расплава даны в работах Cat U4ik&tS а [з] и Острача 4І.
Движения при росте кристаллов разделяются условно на вынужденные, создаваемые подводом механической энергии извне, естественно-конвективные, связанные с действием движущих сил, возникающих вследствие неоднородности расплава (в поле силы тяжести) и, наконец, движения, вызванные специально создаваемыми полями массовых сил (электромагнитные силы). Важную роль в перемешивании расплава играет естественная конвекция, вызванная неоднородностью расплава и обусловленная неравномерностью нагрева и непрерывным изменением концентрационного состава. При этом наличие свободной поверхности и возникновение на ней градиентов сил поверхностного натяжения на счет тлеющейся неоднородности распределения температур и концентрации примесей, может являться причиной поверхностных механизмов движения (термокапиллярная и капиллярно-концентрационная конвекция).
Особенности выращивания кристаллов методом Чохральского рассмотрены Конаковым и др. [б] , Шашковым [б] , а также &{екяег сш НиЯгг ем Н и М&б#а [8] . На рис. I показана схема выращивания кристаллов кремния методом Чохральского по данным работы %/ eeAttt а , Ни&еъ. а [7І .
Метод Чохральского заключается в вытягивании кристалла с помощью затравки из тигля, заполненного жидким расплавленным металлом; для обеспечения тепловой и осевой симметрии применяется вращение кристалла и тигля (см. рис. і). В зависимости от материала имеются аппаратурные отличия в методе Чохральского. Например, для выращивания кристаллов кремния с температурой кристаллизации 1420°С используют кварцевый тигель с резистив-ным нагревом и вращение тигля для тепловой симметрии %ueA/t&t ,
Hu4et (?J . Выращивание гранатов осуществляется из металлических тиглей (платина, иридий) с индукционным нагревом и благодаря хорошей теплопроводности тигля не требует специального вращения тигля ( Мй Ьш & [в] ). Упомянутое отличие вызвано более высокой температурок кристаллизации гранатов (1730 С) и необходимостью использовать дорогостоящие, но обладающие таким важным свойством металлические тигли. Другие аппаратурные отличия имеются при выращивании методом Чох-ральского разлагающихся в обычных условиях материалов, таких как арсенид галлия. В последнем случае для изоляции расплава от внешней среды на поверхность расплава помещают тонкий слой жидкого флюса, который предотвращает интенсивное разложение в период выращивания (см.,например, книгу Шашкова [б] ). Б ряде случаев совершенство кристаллов повышается при программированном периодическом изменении скорости вращения кристалла в период выращивания (см. работу &(г&2 а;Мії еЬ-/(ки/иАа#ї сі
[9] ), а также при использовании дополнительного внутреннего тигля (см. работу Мс\ Ь,ікО! [Ю] ).
Общей задачей технологии является получение совершенных кристаллов с однородной по объему монокристаллической структурой и заданными свойствами (электропроводность, оптические и акустические свойства). На устойчивость монокристаллического роста оказывает влияние кривизна фронта кристаллизации, являющаяся основной технологической характеристикой при получении,например, таких материалов, как арсенид галлия, гадолиний-галлиевые гранаты. Другой важной характеристикой, определяющей однородность электропроводимости, оптических и акустических свойств, является распределение примесей в выращенных кристаллах. Наибольшее значение последний фактор имеет для полупроводниковых материалов (кремний,германий), использующихся для изготовления больших и сверхбольших интегральных схем (БИС, СБИС), а также кри сталлов гранатов для лазерной техники, оптические свойства которых зависят от уровня содержания и распределения примесей.
Наряду с разработкой технологии получения новых материалов в технологической практике существует тенденция к усовершенствованию существующих в настоящее время установок для увеличения объема и габаритов выращиваемых кристаллов. При этом увеличиваются диаметр и длина монокристаллических слитков. Наибольшее продвижение в этом направлении достигнуто при выращивании крупных монокристаллов кремния, диаметром 150-200 мм и длиной 1-І.5 м. Важное значение при этом приобретают результаты математического и лабораторного моделирования, использование которых, наряду с опытом и интуицией специалистов, занимающихся выращиванием кристаллов, позволяет наметить рациональные пути для усовершенствования установок выращивания и выбора технологических параметров.
Обзоры зарубежных работ,выполненных к 1981 г. по моделированию гидродинамических процессов для метода Чохральского, даны hatbdbLs cut [її] и Ниъве. [12] . В этих работах обсуждаются методические вопросы моделирования, а также приводятся данные параметрических исследований отдельного и совместного влияния различных механизмов движения (вращение кристалла и тигля, тепловая, термокапиллярная конвекция и т.д.).
К числу основных направлений численных исследований следует отнести поиск эффектов температурного и концентрационного расслоения и выявление благоприятных факторов для роста кристаллов (см. работы Дубовик, Никитин, Полежаев, Федюшкин [із] и Lanq-OiS (l4j ); а также разработку специализированного математического обеспечения для программированного ведения процесса выращивания ( ІСІУН. и др. [15] ). Важным направлением параметрических исследований является также изу чение влияния изотермических течений расплава, вызванных вращением кристалла и тигля, позволяющее перейти от простых форм движений к более сложным, возникающим в неизотермическом расплаве. Отметим, что этот предельный случай имеет также самостоятельное практическое значение и в условиях невесомости, где гравитационная составляющая отсутствует. Перспективы и конкретные установки для выращивания кристаллов методом Чохральского в космосе рассмотрены в работе Венцля и др. [іб] .
В диссертации отражено современное состояние работ по моделированию метода Чохральского с учетом советских и зарубежных данных и численно исследуются процессы гидродинамики, тепло- и массообмена в расплаве. Основное внимание уделяется методическим вопросам, связанным с разработкой адекватных математических моделей и численных методов, сопоставлением результатов расчетов с данными лабораторных и технологических экспериментов, а также исследуются основные закономерности изотермического и неизотермического течения расплава, его влияния на теплообмен и распределение примесей в объеме и в особенности в подкристальной области. На основе проведенных исследований предлагаются пути рационального выбора технологических параметров для выращивания совершенных кристаллов.
1.2. Состояние вопроса
1.2.I. Вопросы моделирования
К началу работы автора над диссертацией в 1977 г. имелось небольшое количество публикаций по моделированию гидродинамических процессов для метода Чохральского. Недостатком большинства предшествующих работ являлось использование математических моделей, основанных на записи исходных уравнений в приближении пограничного слоя, имеющих ограниченный диапазон применимости и не учитывающих всей совокупности гидродинамических параметров (см. Конаков и др. [б] ). Наибольшее продвижение в разработке более полной модели на основе стационарных уравнений Навье-Стокса было достигнуто в работах
Однако расчеты, проведенные в этих работах, были выполнены для небольшого диапазона малых чисел Рейнольдса и Грасгофа. Основными механизмами, учитывавшимися в расчетах, были вынужденное движение, вызываемое вращением кристалла (тигля) и тепловая гравитационная конвекция, которые не соответствовали режимному диапазону параметров выращивания кристаллов. Отметим также, что в рамках этой стационарной модели нельзя было исследовать нестационарные и колебательные режимы течения и теплообмена, наблюдаемые на практике. Упомянутые недостатки были связаны с ограниченными возможностями применявшейся конечно-разностной методики в целом, аналогичной методике Госмена и др. [18] .
Первые шаги в исследовании изотермических течений расплава при противовращении кристалла и тигля в режимном диапазоне параметров на основе численного решения нестационарных уравнений Навье-Стокса были сделаны LQn.(bu?Ls cu/ [19] . Однако это был единичный расчет и других данных по параметрическим исследованиям не было. По лабораторному моделированию были выполнены работы, дававшие лишь качественное представление об изотермическом течении расплава при вращении кристалла, тигля, например, Туровский, Мильвидский [20] , ChlUtihjzbs , Mssau [2l] , Гришин, Ремизов, Казимиров, Федулов 22] , а также работы, содержащие экспериментальные профили скоростей для изотермических течений -Конаков, Третьяков [23] , и количественные данные о положении поверхности, разделяющей потоки двухвихревого течения при совместном действии вращения кристалла и тепловой конвекции ShitOid [24] . Однако выбор водных растворов в качестве модельной жидкости и красящих веществ, плавающих шариков для визуализации, не позволяет считать эти эксперименты в достаточной мере количественно точными из-за значительного влияния капиллярно-концентрационных механизмов на поверхности жидкости.
Недостатком предшествующего этапа работ являлось также отсутствие сопоставлений данных численного и лабораторного моделирования для метода Чохральского. В этом отношении наиболее интересна работа fioh.so [25] , содержащая результаты такого сопоставления для близкой задачи - течение изотермической жидкости при вращении крышки в цилиндрическом сосуде.
Отметим, что в предшествующий период методические вопросы численного исследования течений при действии тепловой конвекции в диапазоне больших чисел Грасгофа разрабатывались Грязно-вым, Полежаевым [26] . Основными элементами методики этих авторов являлись: численное решение нестационарных двумерных уравнений Навье-Стокса в (w}Ч?) переменных на основе аппроксимации по неявной разностной схеме; новый способ аппроксимации граничных условий для Функции вихря на основе подправления функции тока в приграничных узлах; монотонная аппроксимация конвективных членов по схеме Самарского [27] и итерационное решение уравнения для функции тока с использованием оптимальных итерационных параметров, определенных по Васшпрессу (см. ссылку в работе Шедоренко [28] ). Несколько позже эта методика была усовершенствована для расчета турбулентных режимов конвекции в работе Дайковского, Полежаева, Федосеева [29] . В упомянутых работах также проведены тесты численных решений и сопоставления с данными лабораторного моделирования тепловой конвекции в вертикальных слоях и показана эффективность численной методики при расчете нестационарных режимов в широком диапазоне чисел Грасгофа. Другие разностные схемы и методы решения рассмотрены также в книге Роуча [30] . Рассмотрим работы, появившиеся в период 1977-1984 гг, соответствующий работе автора над диссертацией.
В последние годы методы численного решения значительно усовершенствованы, в связи с чем появились возможности параметрических исследований, превосходящие в отдельных случаях возможности физического моделирования (см. книгу Пасконова, Полежаева, Чудова [31] ). Поэтому все большее значение приобретает проверка адекватности математической модели реальной физической модели. Такая проверка в особенности важна потому, что в реальных условиях при больших числах Рейнольдса и Грасгофа течение теряет симметрию, становится неустойчивым; при этом могут существенную роль приобретать трехмерные эффекты, достоверное описание которых в настоящее время еще затруднительно. Для построения количественных моделей процессов гидродинамики, тепло- и массообмена и распространения их на реальный диапазон режимных технологических параметров важное значение имеют достоверные данные по основным теплофизическим свойствам питающей среды. Обзоры данных, относящихся к расплавам полупроводниковых кристаллов, содержатся, в частности, в работе Глазова, Регеля [32] .
В упомянутый период число работ по моделированию гидродинамических процессов для метода Чохральского резко возросло. Наряду с работами по усовершенствованию метода конечных разностей (МИР) разрабатывался метод конечных элементов (МКЭ), позволяющий учесть геометрические особенности течения в сложных областях с фазовым переходом на границе кристалл-расплав, а также разрабатывались полуаналитические методы решения уравнений Навье-Стокса, позволяющие получить точное численное решение для течений вращающейся жидкости.
Метод конечных разностей (МКР). В работе Ремизова, Смирнова [33] рассмотрено влияние течения изотермического расплава при малых числах Рейнольдса 20, вызванное вращением кристалла и тигля, на распределение легирующей примеси при условии её подпитки на стенках тигля. Ограничение по числу Рейнольдса было вызвано использованием этими авторами методики Госмена и др. [18] и связанными с ней вычислительными трудностями, не позволившими авторам достигнуть режимного диапазона параметров. В работех Люмкиеа, Мартузане, Мартузана [34] , Старшиновой, Фрязинова [35] рассматриваются нестационарные постановки задач, что,в принципе,открывает возможности перехода к диапазону режимных параметров, соответствующих росту объемных монокристаллов большого диаметра. В работеМіск&Ноіс а и др. [36] рассмотрены нестационарные режимы с переменной во времени скоростью вращения тигля. При расчетах распределения примеси в работах Ремизова, Смирнова [33] , Старшиновой, Фрязинова [35] уравнение конвективной диффузии рассчитывается отдельно от уравнений движения, что возможно в пренебрежении концентрационной конвекцией. Построение консервативных конечно-разностных схем для данного класса задач рассматривается в работах Фрязинова [38 и LanfioLs ci [39] . Численное исследование течений при числах Re - 10 и G t 10 связано с преодолением значительных методических и вычислительных трудностей в связи с гидродинамической неустойчивостью и переходом к турбулентности. Только в недавнее время, причем лишь ДЛЯ моделей более общего назначения получены нестационарные численные решения, имеющие стохастические характеристики (Дайковский, Полежаев, Федосеев [29] - тепловая конвекция в вертикальном слое, Шуман, Грецбах, Кляйзер [40] - тепловая конвекция в горизонтальном слое, задача Рэлея-Бенара). Методика и комплекс
программ,предназначенные для такого рода исследований, описаны в работе Бунэ, Грязнова, Дубовика, Полежаева [41] (см.также Пасконов, Полежаев, Чудов [Зі] ).
Методика численного решения уравнений Навье-Стокса на нерегулярных конечно-разностных сетках предложена Фрязиновым J42J ; совместное решение уравнений Навье-Стокса и уравнений фазового перехода (задача Стефана) выполнено в работе Бакировой, Фрязино-ва [43] для вертикального метода Бриджмена. Учет гидродинамики в подкристальном столбике расплава для метода Чохральского выполнен Старшиновой, Фрязиновым [35] . Отметим, что метод конечных разностей благодаря универсальности и сравнительной простоте реализации для задач с прямолинейными границами получил наибольшее распространение. Однако, наряду с упомянутыми преимуществами, имеются трудности при распространении этого метода на задачи с границами сложной формы. Поэтому в последнее время все большее внимание уделяется разработке метода конечных элементов, который отличается от метода конечных разностей специальным выбором аппроксимирующих функций,вариационной формулировкой задачи и возможностью проводить расчеты для сложных областей течения.
Метод конечных элементов (МКЭ). Наибольшее распространение МКЭ получил только в последнее время благодаря появлению современных ЭВМ, обладающих достаточно большим быстродействием и оперативной памятью и позволяющих эффективно решать более сложные системы алгебраических уравнений. Последнее является особенно важным ввиду того, что решение задач методом конечных элементов сводится к решению алгебраической системы уравнений с разреженной матрицей ленточного типа, что значительно усложняет процесс решения по сравнению с МНР, где матрица имеет трех-диагональную структуру, эффективно разрешаемую формулами про гонки. Вопросы численного решения больших разреженных систем уравнений рассматриваются в книге Джоржа, Лю [44] , где основное внимание уделяется предварительной подготовке структуры и оптимизации ширины ленты матрицы.
Теоретические основы использования метода конечных элементов в механике жидкости изложены в книге Коннор, Бреббиа [45] , где рассматриваются вопросы аппроксимации функций на конечных элементах, предлагаются различные вариационные формулировки уравнений МКЭ, а также анализируется точность численных решений и приводятся результаты решения этим методом задач гидродинамики. Консервативные формулировки уравнений МКЭ рассматривались в работах Jesp tSeft- а [46] , L&e., Gxesfu? yCkan а }а/гі [47] . Совместное решение уравнений Навье-Стокса и уравнений фазового перехода (задача Стефана), а также исследование их влияния на распределение примесей выполнено в работе Ctu K.Q Q } fetowh. Q [48] для вертикального метода Бриджмена. Отметим, что в этой работе на основе параметрических исследований МКЭ значительно продвинуто исследование вопросов концентрационного расслоения и радиальной неоднородности распределения примеси на фронте кристаллизации.
В работе Полежаева, Федосеева [49] разработана методика численного решения задач гидродинамики, тепло- и массообмена по МКЭ. Основными элементами этой методики являются: аппроксимация компонент скорости квадратичными и давления - линейными базисными функциями; вариационная формулировка уравнений МКЭ по методу Галеркина и решение алгебраической системы уравнений фронтальным методом (см. работу Hoocf Q [50] ). Большое внимание в работе Полежаева, Федосеева [49] уделено разработке программ для генерации сетки конечных элементов и вопросам оптимизации ленточной структуры матрицы для системы алгебраических урав нений. Кроме этого, в этой работе содержатся данные сопоставления результатов, полученных по МКЭ и МКР, для некоторых задач гидродинамики и делается вывод об эффективности метода конечных элементов.
Применительно к методу Чохральского в работе игооКеі О W&ut& cs c( и др.[5і] по МКЭ проведены расчеты изотермического и неизотермического течения в случае прямолинейных границ области течения. В этой работе рассмотрено течение при отдельном вращении кристалла, противовращении кристалла и тигля, а также при отдельном и совместном действии тепловой конвекции и вращении кристалла в диапазоне изменения числа Рейнольдса до 3 7 10 , числа Грасгофа до 10 при малом значении числа Прандтля -2 равным 10 . Основное внимание эти авторы уделили сопоставлению результатов расчетов с известными численными решениями, полученными по МКР, и установили удовлетворительное совпадение с этими данными для стационарных течений. Кроме этого, при рас А У четах тепловой конвекции для чисел Грасгофа - 2x10 , 10 этими авторами обнаружены колебательные режимы течения, связанные с неустойчивостью течения в пограничном слое вблизи боковой стенки тигля. Анализ температурных полей в кристалле на основе численного решения по МКЭ уравнения теплопроводности совместно с уравнениями фазового перехода выполнен в работе WiMlQIUS Q , HeLcSSet О [52] . Отметим, что на данном этапе возможности моделирования метода Чохральского на основе МКЭ используются не полностью. Например, в упомянутых работах отсутствуют данные расчетов в областях более сложной геометрии, имеющие важное практическое значение (двойной тигель, формообразователь в расплаве).
Другие методы. Наибольшее значение из небольшого числа аналитических решений уравнений Навье-Стокса для моделирования течения в модели метода Чохральского имеет решение для течения, вызванного вращением бесконечного диска в неограниченной жидкости, приведенное в книге Бетчеллора [53] . Это решение на предшествующем этапе использовалось для анализа процессов теплообмена в пограничном слое вблизи вращающегося кристалла (см. книгу Конакова, Веревочкина, Горяинова и др. [5j ). Б последнее время в связи с появлением более адекватных математических моделей метода Чохральского на основе численного решения полных уравнений Навье-Стокса аналитические решения стали использоваться для тестов численных решений, получаемых по МКР и МКЭ. Например, в работах МеМ.ог Q , Chapp Z , btOkQg a [54] , Козельской, Люмкиса [55] получены аналитические решения для течения жидкости между бесконечными вращающимися дисками. При этом точное решение в случае дисков, вращающихся в противоположных направлениях с одинаковыми скоростями, использовано Козельской, Люмкисом [55] для тестов различных аппроксимаций функции вихря на границе области.
Отметим, что важную роль при теоретическом изучении процессов гидродинамики и тепломассообмена при росте кристаллов могли бы играть методы теории устойчивости, позволяющие проводить более полное параметрическое исследование. Однако они, как правило, применяются лишь в простейших случаях при наличии той или иной периодичности или симметрии течения (см.,например, Ckandxasekkat [56] , Гершуни, Жуховицкий [57] ). Некоторые частные случаи течения вращающейся жидкости в замкнутой области рассмотрены в книге Гольдштика [бО] . Таким образом, для численного исследования течения в модели метода Чохральского являются особенно важными данные лабораторных экспериментов по устойчивости течений, являющиеся в настоящее время единственным критерием адекватности математической модели.
Лабораторное моделирование. В последнее время по лабораторному исследованию гидродинамических процессов для метода Чох- ральского также появился ряд новых работ. Наряду с качественным изучением структуры течения при действии различных механизмов движения, в том числе таких, как термокапиллярная конвекция, ранее не учитывавшейся при анализе экспериментальных данных (см. QhitOJCL [24] ), в этих работах начато систематическое количественное изучение характеристик течения и теплообмена. Вопросы лабораторного моделирования применительно к методу Чохральского и конкретные результаты измерений и визуализации течения с помощью модельных жидкостей рассмотрены в работах Уопев я [59] , а/х/эхгок я , Sok&ct&e , Ьскаъжопл а , Seku&Uiss а [60] , JYteoeoo a , IiiW a , Pasksu/ я [61,62] , Бердникова, Борисова
[63] и Бердникова, Борисова, Данченко [б4] .
Добиться равенства критериев подобия в лабораторной модели и в реальной технологической установке даже с точностью до порядка крайне трудно, с одной стороны, в связи с тем, что параметры, которые могли бы изменяться в моделирующей установке, входят в критерии подобия в различной степени, с другой стороны, в связи с сильным отличием в свойствах веществ.Например, использовавшаяся в ряде работ замена расплава германия (вязкость
V 1.3x10 с&г/сек, число Прандтля Pi = 0.016, число Шмидта с?=Ю) водой комнатной температуры ( - I0" WyceK, Pt = 7, SQ =1000) приводит к существенному нарушению подобия тепловых и концентрационных полей в виде различий в числах Н% и So , Если не принимать во внимание это обстоятельство и пытаться реализовать лишь динамическое подобие (подобие полей скорости), то для сохранения равенства чисел Рейнольдса требуется увеличить характерный размер модели, пропорционально корню из вязкости, но при этом может нарушится равенство чисел Грасгофа и т.д. В настоящее время распространенным методом экспериментального исследования стало так называемое "частичное " моделирование,наиболее продвинутое в работах Бердникова, Борисова [63] , Бердникова, Борисова, Панченко [64] для изотермических течений, вызванных отдельным вращением кристалла, в которых также начато изучение неизотермических течений при совместном действии вращения кристалла и тепловой гравитационно-капиллярной конвекции.
Из близких работ, имеющих значение для тестов численных решений, отметим лабораторные и численные исследования изотермического течения жидкости, заключенной в цилиндрический сосуд и приводимой в движение противовращением этого сосуда и крышки на поверхности жидкости ftljkSttQ % HeLj$t а [бб] , а также количественные данные, содержащиеся в работе аУг-Ш1наз о9 FOW&LS a , PtaoSGJc Q , Лее [бб] , где изучены нестационарные течения жидкости, вызванные ускоренным вращением цилиндрического сосуда ( Spin- - Lip ). Применительно к методу Чох-ральского теоретические оценки нестационарного течения изотермической жидкости для режима Sp ii - up сделаны в работе
Математические модели переноса примеси
Рост числа публикаций по численному моделированию гидродинамических процессов для метода Чохральского в 1977-1984 гг. сопровождался появлением новых конечно-разностных методик на основе нестационарных уравнений Навье-Стокса, позволяющих проводить исследования в режимном диапазоне параметров выращивания кристаллов; разработкой метода конечных элементов, позволяющего проводить исследования в областях со сложными границами, а также разработкой и использованием для тестов приближенных методов, позволяющих получить точное решение. Расширение возможностей численного моделирования позволило перейти к расчетам термокапиллярных эффектов, структуры течения при ускоренном и замедленном вращении кристалла и тигля, а также параметрическим исследованиям в режимном диапазоне параметров выращивания кристаллов и сопоставлениям ( появившимися в последнее время данными по лабораторному моделированию метода Чохральского.
Однако в упомянутых работах вплоть до настоящего времени остаются нерешенными вопросы, связанные с количественной проверкой адекватности численных решений по данным лабораторного моделирования; вопросы устойчивости течений, а также ввиду многопараметрической обусловленности задачи остается неисследованным влияние таких важных факторов, как температурное профилирование на стенках тигля, изменение скоростных режимов врашения кристалла и тигля, зависимость течения и теплообмена от отношения радиуса кристалла и тигля, вопросы, связанные с влиянием течения на распределение примесей (кислород, легирующая примесь) и т.д.
Весь период работы автора над диссертацией можно разделить на два следующих этапа: 1977-1980 гг. и І98І-І984 гг. На первом этапе (1977-1980 гг.) была разработана математическая модель гидродинамических процессов для метода Чохральского и численная методика на основе МКР, выполнены тесты численных решений путем сопоставления с данными, полученными по другим разностным схемам, и данными лабораторного эксперимента W&1H- -УаЧ.па$а и др. [бб] . В это же время исследовано влияние отдельных механизмов движения и начато параметрическое изучение изотермических и неизотермических течений. На втором этапе (І98І-І984 гг.) была усовершенствована конечно-разностная методика: осуществлен переход к расчетам с двойной точностью, уменьшены затраты машинного времени и повышены аппроксимационные свойства решения (аппроксимация вихря на границе полной области, введение релаксации по граничным условиям для вихря); разработана математическая модель метода Чохральского на основе МКЭ; проведены тесты с аналитическими и численными решениями; выполнены сопоставления с данньши лабораторного эксперимента; проведены параметрические исследования изотермических и неизотермических течений в широком диапазоне параметров; исследовано влияние течения расплава на теплообмен и распределение примесей; предложены пути рационального выбора технологических параметров, которые были апробированы в заводских условиях и дали положительные результаты.
Результаты диссертационной работы докладывались на 38-й, 39-й конференциях Латв. гос. университета (Рига, 1978 и 1979); УШ и IX Всесоюзных семинарах по численным методам механики вязкой жидкости (Томск, 1979; Ленинград, 1982); семинаре Института теплофизики СО АН СССР (Новосибирск, 1979); УІ Международной конференции по росту кристаллов (Москва, 1980); Ш и ІУ Университетских школах "Нелинейные задачи гидродинамической устойчивости" (Ко любакино, 1980 и 1982); І-Ш Всесоюзных семинарах по гидромеханике и тепломассообмену в невесомости (Москва,1979; Пермь, 1981; Черноголовка, 1984); У Всесоюзном съезде по теоретической и прикладной механике (Алма-Ата, 1981); Научно-технических семинарах в Гиредмете (Москва, 1980, 1981, 1983); семинаре Института кристаллографии АН СССР (Москва, 1981); УП Всесоюзной школе "Теоретические и прикладные проблемы вычислительной математики и математической физики" (Рига, 1982); 1-м Всесоюзном семинаре "Математическое моделирование процессов затвердевания металлов и сплавов" (Новосибирск, 1983); семинаре ИТР завода ОХМЗ (Подольск, 1984); семинарах Института проблем механики АН СССР (Москва, 1984).
Основные результаты диссертации опубликованы в 16 работах [96 - ІІІІ и применялись в Институте проблем механики АН СССР при проведении плановых работ по комплексной проблеме ГКНТ СССР (шифр 0.09.II), договору о социалистическом содружестве с ИТФ СО АН СССР. Результаты диссертации используются в
ГИЕЕДМЕТё и ОХМЗ МЦМ СССР, НИИМВ МЭД СССР, ИТФ СО АН СССР при проведении технологических и лабораторных экспериментов с целью выбора оптимальных параметров выращивания кристаллов. Комплекс программ по методу конечных разностей внедрен в ИФИ АН Арм, ССР и используется при параметрическом исследовании технологических режимов. В приложении приведены документы о внедрении результатов диссертационной работы в ИФИ АН Арм. ССР и на заводе ОХМЗ МЩ СССР. Автор выражает глубокую благодарность своему научному руководителю В.И. Полежаеву за постоянное внимание и помощь в работе.
Метод конечных элементов
Математическое моделирование метода Чохральского основано на численном решении полных, нестационарных уравнений Навье-Стокса совместно с уравнениями переноса тепла и примеси в приближении Буссинеска и предположении осевой симметрии. Математическая модель учитывает основные механизмы движения: вращение кристалла и тигля со скоростями П% и fLy t соответственно, гравитационную тепловую и концентрационную конвекцию, а также поверхностные механизмы движения - термокапиллярную и капиллярно-концентрационную конвекцию. Понижение уровня расплава и увеличение радиуса кристалла в период выращивания учитывается путём проведения расчетов для соответствующих значений геометрических параметров H/RT , RK/@T (см.рис. 2- и 3). Кроме этого, изменение характера теплообмена учитывается путем перераспределения потоков тепла или профиля температур на стенках тигля, свободной поверхности и кристалле. Полная классификация механизмов движения расплава приведена на рис.4.
Таким образом, в математическую модель включено много различных факторов, влияние которых на течение пока недостаточно изучено (см. обзор литературы во введении в п. 1.2). Поэтому на данном этапе при моделировании не учитывалось тепловыделение при кристаллизации (задача Стефана) и предполагались фиксированными граница кристалл-расплав и поверхность жидкости. Отметим, что такой подход является приемлемым при условии, что температурные граничные условия для численного моделирования задаются по экспериментальным данным.
В диссертации для модели метода Чохральского численное исследование процессов гидродинамики и тепломассообмена в одно-связных областях с прямолинейными границами (рис. 2) осуществляется методом конечных разностей и в областях сложной формы -(рис. 3) - методом конечных элементов. В зависимости от применяемого численного метода для формулировки исходных уравнений и граничных условий использованы следующие переменные: (w г , W } & ) - для метода конечных разностей, ( U , , /Ъ } W f & ) -для метода конечных элементов. При этом отличия в выборе переменных объясняются стремлением использовать при решении уравнений Навье-Стокса наиболее эффективные аппроксимации граничных условий в рамках каждого из упомянутых методов. Эти переменные приводятся к безразмерному виду;в качестве масштабов взяты следующие величины: радиус кристалла R , окружная скорость вращения кромки кристалла ліц Г\ц » масштаб времени 4/11« , максимальный перепад температур между стенками тигля и кристаллом д 7-7 = Тт - TL и разность концентраций примеси на стенках тигля и кристалле Л = Ст Сц , Полная сводка параметров математической модели и формул для их вычисления приводится в списке "Условных обозначений" в приложении.
Схема математической модели для численного исследования на основе метода конечных разностей показана на рис. 2 и предполагает односвязность области течения и прямолинейность её границ. Граничные условия соответствуют прилипанию жидкости на стенках тигля и поверхности кристалла, скольжению жидкости вдоль свободной поверхности с учетом термокапиллярных эффектов, а также заданным условиям теплообмена на всех участках границы. Таким образом, граничные условия записываются в следующем виде
Схема математической модели на основе метода конечных элементов показана на рис. 3 и предполагает возможность сложных границ области течения: сферическое скругление дна тигля, криво-линейность границы расплав-кристалл, двойной тигель и т.д. С учетом осевой сишлетрии исходная нестационарная система уравнений Навье-Стокса совместно с уравнением переноса тепла в - переменных и в приближении Бус-синеска записывается в следующем виде
Математические модели переноса примесей в настоящее время разработаны в меньшей степени, чем модель конвективного теплообмена. Последнее связано с отсутствием аналогичных данных по распределению концентрации примесей на границе расплава. Поэтому в настоящее время математическое моделирование процесса переноса примесей позволяет сделать лишь оценку влияния течения на содержание и однородность распределения примесей на фронте кристаллизации. При этом для интерпретации результатов расчетов важную роль играют контрольные измерения распределения примесей в кристалле, соответствующие фиксированным условиям процесса выращивания этого кристалла, что позволяет установить диапазон адекватности математической модели. Наибольшие продвижения в этом направлении сделаны при исследовании распределения легирующей примеси и кислорода для расплава такого простейшего полупроводникового материала, как кремний.
Математическое моделирование переноса примеси выполнено на основе совместного численного решения гидродинамической части задачи, сформулированной в предыдущем параграфе 2.1, и нестационарного уравнения конвективной диффузии, записываемого в следующем виде
При.этом при исследовании нестационарных режимов переноса примеси и оценке совместного действия тепловой и концентрационной конвекции осуществлялось одновременное решение уравнений гидродинамики и переноса примеси.
Противовращение кристалла и тигля
При изменении / от 1240 до 3300 последовательность мгновенных картин линий тока, соответствующих перестройке течения через безразмерный интервал времени, равный A t =50 от начального времени t0 =100, показана на рис.24 а-е. В этом случае отчетливо наблюдается влияние перестройки, происходящей в основном потоке, на конфигурацию линий тока и положение вторичного вихря в подкристальной области. В начальный момент =100 (рис. 24 а) вторичный вихрь расположен вблизи кристалла, в следующий = 150 (рис. 24 б) происходит деформация основного вихря,и он смещается ко дну тигля, затем при Ч1 200, 250 (рис. 24 в,г) в связи с колебательным характером изменения структуры основного потока вторичны;", вихрь опять поднимается вверх, и, наконец, при 2 5=350, 450 устанавливается стационарная картина течения с положением вторичного вихря на расстоянии 1/ЗхК от дна тигля (рис. 24д,еУДля диапазона больших чисел Рейнольдса при переходе от Re= 3300 к 1.57x10 перестройка течения показана на рис. 25: (а) мгновенная картина линий тока в начальный момент времени, характеризующаяся усилением радиального течения вблизи свободной поверхности и смещением центра основного вихря в угловую область; (б) промежуточная картина течения, соответствующая перестройке основного потока и образованию вторичного вихря вблизи боковой стенки тигля; и (в) окончательная стационарная структура течения с вторичными вихрями вблизи оси и боковоіі стенки тигля (±тах 2.8x10 ).
Для обобщения всех режимов течения на рис. 26 приведена диаграмма в плоскости ( H/Rj t Re ), разделяющая различные типы течений (линии 1,П,Ш и 1а, Па, Ша): появление и развитие основного течения, возникновение вторичного вихря, колебательные режимы течения. Эта диаграмма впервые построена в работе Бердникова, Борисова [бЗ] на основе экспериментальных данных и в диапазоне существования стационарных течений удовлетворительно согласуется с результатами численных расчетов, проведенных для случая плоской формы диска (кристалла) без наличия капиллярного поднятия жидкости. На диаграмме звездочкаїли () показаны точки, соответствующие проведенной серии численных параметрических расчетов. На диаграмме в координатах: число Рей-нольдса - Re и относительная высота слоя жидкости в сосуде (тигле) - H/RT » систематизированы результаты исследований пространственной формы течения при фиксированном отношении Rc/T= о.Зб.
Точки 1-3 соответствуют модели с плоским диском (кристаллом); точки 4-6 соответствуют модели с выпуклым диском и точки 7-9 - модели с вогнутым диском. Дополнительные точки 10,11 соответствуют модели с плоским диском при наличии капиллярного поднятия жидкости на максимальную высоту п./RT = 0.11. В последнем случае наличие подкристального столбика жидкости приводит к затягиванию развития течения, и границы На и Ша соответствуют такой же пространственной организации течения, что и границы П и Ш, но без капиллярного столбика жидкости. Кроме этого, как видно из рис. 26, кривизна диска не оказывает существенного влияния на границу развития течения по числу Re и не искажает существенно линии тока.
По точкам слева от кривой I течение имеет двухярусную структуру ( H/Rj , Re ). Вблизи оси жидкость поднимается к диску, затем отбрасывается к боковой стенке тигля, опускается и вновь подтягивается вблизи оси в виде подъемной струи. При
Ке R&i Б описанном течении участвует только часть жидкости, находящейся вблизи поверхности диска. При этом в эксперименте в придонной области наблюдается вторичное, очень слабое движение жидкости в обратном по компонентам скорости щ iF направлении. Однако с ростом Re основное течение занимает все больший объем и при нижняя область с обратным течением исчезает.
В области от кривой I до кривой П структура течения остается одновихревой, но пространственная ориентация основного вихря изменяется (см. п. 4.1.I). С ростом Re растет скорость циркуляции жидкости. На границе П в подкристальной области появляется приосевая особенность, которая с ростом числа Re развивается в своеобразное вторичное течение, форма которого зависит от геометрического параметра Н/г?т (см. п. 4.1.2).
На границе в эксперименте течение приобретает колебательный- характер. Слева от границы течение стационарное и осесимметричное. Вид колебаний при Re Re l зависит от параметра Поиски аналогичной границы расчетным способом осуществлялись на основе ІЖР до значительно больших, чем в эксперименте значений числа Рейнольдса Re =2.0x10 при максимальном заполнении сосуда, согласно диаграмме Н/т= 1.4. Эти расчеты были проведены на подробной неравномерной сетке с количеством узлов 61x81, сгущенных в областях наибольших градиентов искомых функций. Однако колебательные процессы, которые возникали за счет изменения начальных данных при переходе к большему числу Рейнольдса,со Бременем затухали и решение выходило на стационарный режим. Возможные причины этого явления рассмотрены выше в п. 4.1.3.
Некоторые эффекты совместного действия тепловой конвекции и вращения кристалла
Участок П соответствует периоду выращивания кристалла максимального фиксированного диаметра с развитым двухвихревым течением в тигле, сохраняющим свою структуру и положение "точки встречи" X =0:3 при изменении значения параметра irt- /Re что связано с преимущественно количественным перераспределением интенсивности потоков для достаточно широкого диапазона изменения этого параметра. Отметим, что этот эффект характерный также для изотермического двухвихревого течения при противовращении кристалла и тигля, обнаружен / при Re = 688 и изменении в широком диапазоне параметра Участок Ш соответствует преобладанию течения от вращения кристалла и значительному смещению "точки встречи" к боковой стенке тигля. Последнее может характеризовать конечный период выращивания кристалла, когда происходит значительное понижение уровня, расплава и уменьшается вклад тепловой конвекции.
Отметим, что изменение числа Прандтля не столь существенно отражается на положении "точки встречи" потоков (см. данные, отмеченные на рис. 46,а соответствующими метками) по сравнению с перераспределением распределения температур на стенках тигля. Влияние последнего фактора обнаруживается по поведению пунктирной линии, построенной по данным работы и равномерном распределении температур на стенках тигля. При этом отличия от участка П-Ш сплошной линии связаны с образованием вблизи оси вторичной циркуляции, вызванной донным нагревом тигля и приводящей к ослаблению подкристального течения расплава от вращения кристалла. Влияние термокапиллярных эффектов на "точку встречи", потоков можно оценить по зависимости X =#uV v , приведенной на рис. 46,6 и соответствующей Qt = 8.2x10 , Pt = 15, равномерному нагреву боковой стенки тигля. Анализ этой зависимости показывает, что влияние термокапиллярных эффектов существенно и при увеличении числа Марангони до Л/#, =5.12x10 приводит к значительному перемещению "точки встречи" в сторону кристалла.
Отметим, что установленная зависимость положения "точки встречи" двух потоков X = H tu/Re ) справедлива лишь при стационарной картине течения, что соответствует не слишком большим числам Рейнольдса, Грасгофа и т.д. Проведенные в последнее время лабораторные исследования неизотермических течений для модели метода Чохральского обнаружили переход к регулярным колебаниям температуры при вращении кристалла и действии тепловой конвекции (см. работы fksko- , M/to&tT fet ed \б\ , Бердников, Борисов, Панченко, Простомолотов [ill] ).
Противовращение кристалла и тигля применяется при выращивании большинства полупроводниковых кристаллов и, в частности, такого важного материала, как кремний (см. Шашков [б] ). Постоянное совершенствование технологии выращивания кремния в направлении увеличения диаметра и длины кристаллов сопровождается одновременным повышением требований к однородности выращиваемых слитков, зависящей от распределения легирующей примеси и кислорода. Поэтому основной задачей численного моделирования является поиск таких режимов течения расплава, которые приводят к однородному распределению примесей на фронте кристаллизации. Последнее сводится к параметрическому исследованию функционального соотношения (2.54) и рассматривается далее по тексту для наиболее характерных технологических режимов выращивания кристаллов.
Рассмотрим некоторые особенности неизотермического течения при противовращении кристалла и тигля в случае равномерного и преимущественно бокового нагрева стенок тигля, а также следующего диапазона изменения режимных параметров
Для определенности основным будем считать течение, соответствующее равномерному распределению температур на стенках тигля и /?е =З.ЗхЮ3, Qt =107, У& 1,/%=0, ////. =0.66. Напомним, что отдельное действие Тепловой конвекции в этом диапазоне параметров рассмотрено ранее в п. 5.1.2 и соответствует структуре течения и изотерм, приведенным на рис. 38,а, причем изотермическое течение при противовращении кристалла и тигля рассмотрено в предыдущей главе 4 и показано на рис. 32,а. Сопо - 124 ставление интенсивности течения для этих двух отдельных случаев, а также предварительная оценка ускорений, вызванных вращением кристалла и тигля (В -О-к х =8.83 см/сек и $Т =/2-х у- = 0.33 см/сек ), и их сопоставление с ускорением силы тяжести Q, = 980 см/сек: , позволяет сделать вывод о преобладании тепловой конвекции в тигле при их совместном действии. При этом заметное влияние вынужденной конвекции возможно лишь вблизи вращающейся поверхности кристалла. Рассмотрим результаты численных расчетов для совместного действия тепловой конвекции и противовращения кристалла и тигля. На рис. 47 показана структура изотерм и линии тока при Ми =0 - (а) и Мп. =10 -(б), соответствующие в обоих случаях двухвихревому течению расплава. При этом наиболее интенсивное течение возникает вблизи боковой стенки тигля за счет действия тепловой конвекции. Аналогично случаю отдельного действия тепловой конвекции слои расплава, нагретые вблизи боковой стенки до максимальной температуры, поднимаются вверх вдоль её поверхности, охлаждаются за счет теплообмена с внешним пространством на поверхности расплава и затем опускаются ко дну тигля.
