Содержание к диссертации
Введение
1 Алгоритмы численного исследования фазовых переходов в многокомпонентных системах 18
1.1 Математическая модель 18
1.1.1 Основные уравнения 19
1.1.2 Условия на фронте кристаллизации 20
1.1.3 Начальные данные 21
1.2 Алгоритм решения 21
1.2.1 Замена переменных 21
1.2.2 Разностная задача 25
1.2.3 Метод решения сеточных уравнений 28
1.2.4 Результаты тестовых расчетов. Сравнение с полунеявным алгоритмом 30
2 Моделирование диффузионных режимов выращивания эпитаксиальных слоев . 47
2.1 Изотермический рост 47
2.2 Рост с предварительным подрастворением подложки. CdTe подложка 56
2.3 Влияние диффузии в твердой фазе на динамику процесса ЖФЭ 62
2.4 Рост из недосыщенного раствора. CdxHgi-xTe подложка . 66
3 Математическое моделирование ЖФЭ четверных соединений в двумерном приближении. 76
3.1 Математическая модель 76
3.1.1 Дифференциальная задача 77
3.1.2 Граничные условия на фронте кристаллизации . 78
3.1.3 Начальные данные 80
3.2 Разностная задача 81
3.2.1 Модельный пример 81
3.2.2 Решение концентрационной задачи 84
3.2.3 Решение системы сеточных уравнений. Метод Ньютона 86
3.2.4 Решение линеаризованной задачи 88
3.2.5 Примеры тестовых расчетов 89
3.3 Результаты расчетов 102
Приложение 106
Список литературы 107
- Условия на фронте кристаллизации
- Рост с предварительным подрастворением подложки. CdTe подложка
- Рост из недосыщенного раствора. CdxHgi-xTe подложка
- Граничные условия на фронте кристаллизации
Введение к работе
Математическое моделирование процессов кристаллизации многокомпонентных систем находится в центре внимания исследователей. В значительной степени работы в этой области стимулируются проблемами, возникающими в технологии получения новых полупроводниковых материалов. Высокие требования, предъявляемые современным приборостроением к свойствам и качеству используемых полупроводниковых структур и материалов, могут быть удовлетворены лишь на основе детального изучения физических процессов, происходящих при их получении, и тщательной отработки самого технологического процесса. Определяющую роль в формировании свойств получаемых материалов играют фазовые превращения и процессы тепломассопереноса, формирующие условия на границе раздела фаз. Натурные эксперименты в этой области сложны и трудоёмки, требуют больших энергетических затрат, в них используются дорогостоящие, зачастую ядовитые материалы, а длительность одного эксперимента может достигать нескольких недель. Поэтому использование методов математического моделирования в исследовании и оптимизации технологических режимов и установок, для анализа результатов физических экспериментов является неотъемлемым этапом на пути создания способов получения материалов с заданными свойствами.
Однако, до настоящего времени математическое моделирование в задачах технологии ограничивается в основном рассмотрением процессов кристаллизации в чистых веществах или бинарных растворах-расплавах. Вместе с тем, тепломассоперенос и фазовые превращения в трех и четырехкомпонентных системах лежат в основе многих промышленных способов получения объемных монокристаллов и полупроводниковых структур и представляют значительный теоретический и практический интерес.
Данная работа посвящена численному исследованию процессов кристаллизации в многокомпонентных системах. В качестве примера, иллюстрирующего область применения полученных результатов, рассматривается один из наиболее распространенных способов получения монокристаллических структур - метод жидкофазовой эпитаксии (ЖФЭ), который представляет собой кристаллизацию из раствора-расплава на подложке определенной кристаллографической ориентации. В простейшем случае насыщенный при начальной температуре раствор кристаллизующихся компонентов приводится в контакт с подложкой, и вся система охлаждается по заданному закону. С понижением температуры жидкая фаза становится пересыщенной, и растворенные в ней вещества осаждаются на границе расплав - подложка в виде монокристаллического слоя [1]. В ходе всего процесса на границе раздела фаз имеет место квазиравновесие, т.е. концентрация растворенных компонентов у фронта кристаллизации в любой момент времени близка к равновесному при данной температуре значениям. Перенос растворенного вещества к поверхности раздела фаз осуществляется механизмом диффузии и конвекции [2]. Пересыщение в объёме может создаваться различными способами: с помощью принудительного охлаждения системы подложка -расплав, предварительного охлаждения, то есть до контакта с подложкой, охлаждения жидкой фазы ниже температуры насыщения, подпитки раствора - расплава за счёт растворения подложки, изотермического смешивания равновесных растворов и т.д. [1-3]. Наиболее часто в технологии эпитаксиального выращивания полупроводниковых материалов из жидкой фазы применяется способ, при котором температура системы расплав - кристалл изменяется по заданному закону. Этот вариант ЖФЭ сравнительно прост в техническом отношении, обладает достаточной гибкостью, дает возможность выращивать многослойные эпитакси-альные структуры с контролируемыми параметрами каждого из слоев [1, стр. 36].
В некоторых технологических режимах подложка помещается в ненасыщенный при температуре контакта расплав. Вся система выдерживается некоторое время при постоянной температуре, а затем начинается процесс охлаждения. В этом случае росту эпитаксиального слоя (ЭС) предшествует процесс растворения подложки, благодаря которому происходит обогащение раствора - расплава кристаллизующимся веществом.
Условия на фронте кристаллизации
Предложен новый чисто неявный алгоритм решения одномерной нестационарной задачи кристаллизации многокомпонентных растворов. Учет движения фронта кристаллизации осуществляется с помощью замены переменных типа Ландау [4,21], в результате которой граница раздела фаз остается неподвижной на протяжения всего процесса. Разностная аппроксимация строится в новой системе координат на неподвижной сетке с помощью интегро-интерполяционного метода [64]. В отличие от традиционно используемых подходов, коэффициенты уравнений, связанные с преобразованием системы координат и определяемые скоростью движения границы фазового перехода, берутся с верхнего временного слоя. Дифференциальная задача аппроксимируется с помощью консервативной, чисто неявной схемы. Соответствующая система нелинейных сеточных уравнений решается методом Ньютона относительно вектора неизвестных, компонентами которого являются все искомые концентрации в жидкой и твердой фазе и скорость движения фронта кристаллизации.
По сравнению с известными ранее методами, построенный алгоритм обладает большим запасом устойчивости (является практически безусловно устойчивым) и гарантирует получение надежных результатов для широкого класса задач о фазовых переходах. К достоинству метода надо отнести отсутствие необходимости согласования начальных составов жидкой и твердой фаз на фронте кристаллизации. Совместное определение полей концентраций и нового положения границы раздела фаз решает ключевую проблему выбора начальных данных, существующую при последовательном вычислении распределения состава и положения фронта кристаллизации [4,5,35-42].
На примере моделирования процесса жидкофазовой эпитаксии трех-компонентных соединений проводится сравнение предложенного алгоритма с методом, в котором преобразование сетки осуществляется по значениям искомых функций с нижнего временного слоя. В частности показано, что последовательное определение распределения состава в системе и нового положения фронта кристаллизации может приводить к физически неверным результатам.
Предложенный численный метод решения задач о фазовом переходе в многокомпонентных системах строится на примере ЖФЭ тройных соединений, однако алгоритм легко обобщается на случай кристаллизации растворов с произвольным числом растворенных компонент, а также позволяет учесть выделение или поглощения тепла на границе раздела фаз.
Однородность распределения состава по толщине слоя обеспечивают режимы выращивания, при которых перенос растворенных компонентов к фронту кристаллизации осуществляется преимущественно механизмом диффузии. С точки зрения качества получаемых материалов, такие режимы наиболее предпочтительны и их исследование представляет значительный интерес. Анализ распределения состава в направлении роста предполагает изучение влияния диффузии в твердой фазе. В большинстве работ, посвященных численному исследованию конвективного массопереноса в растворе, диффузия в твердой фазе не учитывается [46-48,55,56]. Такой подход является вполне оправданным, так как этот процесс не оказывает существенное влияние на движение расплава. Однако, диффузия в твердой фазе играет важную роль в формировании профиля концентрации, особенно вблизи гетерограницы. Кроме того, учёт диффузии в твердой фазе позволяет построить самосогласованную математическую модель процесса, описывающую как рост монокристаллического слоя, так и его растворение [58].
Задача решается с помощью процедуры расщепления по физическим процессам [66]: сначала вычисляется поле скоростей, а затем определяется распределение концентраций в системе. В данной главе основное внимание уделяется решению концентрационной задачи. Хорошо известно, что в рассматриваемом классе задач численная реализация граничных условий [13, 31-33,67] определяет устойчивость алгоритма. Были опробованы всевозможные схемы, в которых распределение растворенных компонент определяется последовательно или две из них вычисляются совместно, а затем третья досчитывается. Однако все эти алгоритмы оказались неустойчивым в рассматриваемом диапазоне технологических параметров.
Надежный метод удалось построить лишь на основе совместного вычисления распределения концентраций в системе. Но это сопряжено со значительными трудностями, поскольку на каждом шаге по времени приходится решать систему уравнений относительно вектора, компонентами которого являются сеточные значения концентраций всех растворенных компонент, и свойства соответствующего разностного оператора неизвестно. Это препятствует применению итерационных алгоритмов, требующих для своей реализации априорной информации о спектре. Использование прямых методов решения малоэффективно, так как требует больших затрат машинного времени. Попытки обобщить алгоритмы, применяемые для решения аналогичных задач в случае двух- и трехкомпонентных систем не дали удовлетворительных результатов. Поэтому была разработана новая вычислительная процедура для решения концентрационной задачи.
Рост с предварительным подрастворением подложки. CdTe подложка
Рассмотрим режимы выращивания ЭС в условиях принудительного охлаждения с предварительным подрастворением подложки. Такие режимы широко используются в технологической практике. Растворение подложки перед началом роста благоприятно сказывается на качестве материала, способствуя насыщению раствора, устраняет загрязнения с поверхности подложки. Однако, этот процесс очень сложно контролировать.
Воспользуемся для анализа поведения системы на начальной стадии процесса результатами раздела 2.1 (стр. 47). После контакта подложки с неравновесным раствором система стремится к квазиравновесию, которое может достигаться как за счет растворения подложки, так и за счет роста ЭС. Как именно будет развиваться процесс зависит от внутренних движущих сил ЖФЭ: "пересыщения твердой фазь х — Хо и недосыщения раствора 6Т = Tsat — То. В данном случае они действуют в противоположные стороны. Подложка CdTe является более тугоплавкой, чем солидус исходного раствора Cd Hgi- Te, поэтому "сила" XQO —хо способствует росту эпитаксиального слоя. С другой стороны, раствор недосыщен (Tsat TQ), и значит подложка должна растворяться, насыщая жидкую фазу. Учитывая, что диффузия в твердой фазе мала, а величина ST достаточно большая, то, в соответствии с оценкой (2.1.1), в системе начнется расворение подложки и влияние " пересыщения твердой фазы" (хоо — хо) будет несущественно. Эволюция состава жидкой фазы на гетерогранице. Влияние начального недосыщения раствора 5Т = Tsat То. Кружками обозначены начальные состав раствора для различных 5Т. 1 -изотерма Т = 782.2 ДГ для ST = -2К, 2 - изотерма Т = 781.3 для ST = -ЪК, 3 - изотерма Г = 780.5 для ST = -10К.
Эволюция состава жидкой фазы на фронте кристаллизации в тече-ниее всего процесса изображена на рисунке 2.9. В результате мгновенно установившегося квазиравновесия состав жидкой фазы на границе расплав - кристалл принимает значение, которое лежит на изотерме То, соответствующей начальной температуре системы (рис. 2.9). Необходимо заметить, что равновесие устанавливается за счет изменения состава жидкой и твердой фаз (рис.6 1.16, ST = —10if). Участок кривой от То до Т (рис. 2.9) соответствует периоду растворения подложки, затем начинается рост эпитаксиального слоя (Т Т ). Т - температура, при которой происходит смена растворения ростом. При температуре Т = 11ЪК охлаждение прекращается, и рост продолжается в изотер 6На всех рисунках в твердой фазе приводится распределение только С ґГе, а распределение НдТе легко определяется в соответствии с уравнением (1.1.2) мическом режиме, поэтому состав жидкой и твердой фазы на фронте кристаллизации постоянен [47,59,81].
В рассматриваемом случае подложка не содержит Нд, поэтому при растворении в раствор поступает только Cd (рис. 1.16, 5Т = —10К). В то же время Нд диффундирует в твердую фазу, вблизи фронта кристаллизации содержание Нд в растворе уменьшается. Одновременно с растворением в подложке формируется тонкий диффузионный слой CdxHgi_xTe, который растворяется прежде, чем Нд успевает проникнуть вглубь подложки (рис. 2.10). Концентрация раствора приближается к равновесной. Благодаря растворению подложки и охлаждению системы состав жидкой фазы постепенно выравнивается по пространству (рис. 2.11, 2.12), и под влиянием продолжающегося охлаждения начинается рост ЭС (рис. 2.10).
На развитой стадии процесса распределение состава жидкой фазы монотонно. На фронте кристаллизации содержание Cd и Нд в расплаве Те минимально и определяется уравнением ликвидуса при текущей температуре Т. Таким образом, раствор насыщен на границе раздела фаз и пересыщен в глубине жидкой фазы. Лишь в короткий период времени, когда происходит смена растворения ростом, в расплаве можно наблюдать немонотонное распределение состава (например, линия 4 рис. 2.11). В случае, приведенном на рисунке 2.11, локальное пересыщение невелико и составляет 0.17К. Величина локального пересыщения зависит от коэффициента диффузии в жидкой фазе, толщины расплава, степени недосыщенности раствора 5Т. Сильное локальное пересыщение может повлечь за собой возникновение процесса объемной кристаллизации внутри расплава. В рассматриваемых условиях выращивания увеличение толщины раствора до 2мм приводит к усилению локального пересыщения до 0.35.К\
Для изучения влияния диффузии в твердой фазе на динамику процесса ЖФЭ была проведена серия расчетов: коэффициент диффузии Ds Cd варьировался в диапазоне 5-10-13см2/с- 5-1СГ9см2/с, хо = 0.207, Хоо = 1, 5Т = -2К, Т0 = 783#, а = 0.5#/мин, Dl(cd = D1 = 5 ИГ5см2/с. Рисунок 2.14 демонстрирует распределение состава в системе на ранней стадии процесса для Ds Cd — 5 10_11см2/с, 5Т = —2К. Рисунки 1.16 и 2.14 качественно похожи: в жидкой фазе наблюдается увеличение содержания Cd и уменьшение доли Нд, в твердой - формирование диффузионного слоя CdxHg\-xTe, однако при Ds(Cd) = 5 10_13см2/с (рис. 1.16) происходит растворение подложки, при Ds(CdJ = 5- 10 псм2/с (рис. 2.14 ) - рост эпитаксиального слоя. Направление движения фронта кристаллизации определяется соотношением суммарных потоков Cd и Нд через границу раздела фаз. В случае Ds Cd — 5 10_13см2/с результирующий поток Cd, складывающийся из потоков Cd в твердой и жидкой фазах, также как и суммарный поток Нд, направлен в жидкую фазу, чем обусловлено растворение подложки, а при Ds Cd ) = 5 10_11см2/с -в противоположную сторону, то есть процесс кристаллизации.
Распределение состава в жидкой и твердой фазе на ранних стадиях процесса для различных коэффициентов диффузии в твердой фазе приведено на рисунке 2.15. При сильной диффузии в твердой фазе jjs(Cd) 5 .10-11см2/с) перегрев системы ST == — 2 К уже не определяет ее поведение, и под действием "пересыщения твердой фазы" XQO — х0 начинается рост ЭС.
Рост из недосыщенного раствора. CdxHgi-xTe подложка
В результате устанавливающегося в системе квазиравновесия, состав твердой фазы на границе раствор - кристалл становится равным хе Хо, поэтому внутри подложки образуется тонкий слой, в котором содержание Cd выше, чем х (рис. 2.18 распределение 1). Пока идет процесс растворения структура этого пограничного слоя практически не изменяется, только под влиянием охлаждения уменьшается значение хе на поверхности подложки (рис. 2.18 распределение 2, 3). Рост из недосыщенного раствора на Cdo orHдолэзТе подложку. Распределение Нд в растворе. 6Т = —10К. 1 - время 1с, 2 - время 1мин, 3 - время 5мин, 4 - время 6.5мин, 5 - время 7мин, 6 - время Юмин. фронта кристаллизации раствор становится насыщенным (рис. 2.19, 2.20 распределения 4, 5, 6), начинается рост ЭС. Состав жидкой фазы к этому моменту выравнивается (рис. 2.19, 2.20 распределения 4), но в подложке зона повышенного содержания Cd, возникшая на стадии растворения, остается (рис. 2.18 распределение 4, 5, 6). В результате на гетерогранице подложка - ЭС происходит резкое изменение состава твердой фазы (см. рис.8 1.4). Максимальное значение концентрации CdTe в ЭС определяется температурой начала роста (Т ) и величиной диффузии в твердой фазе.
Распределение CdTe в твердой фазе для различных 8Т приведено на рисунке9 2.21. Величина скачка на гетерогранице в данном случае определяется величиной первоначального переохлаждения раствора. При небольших значениях 6Т равновесный состав хе мало отличается от исходного состава подложки Хоо = хо и "всплеск" концентрации CdTe на границе подложка - слой невелик. Для 5Т = — 2К хтах — Хоо = 0.0003, где хтах -максимальная концентрация CdTe в твердой фазе. С ростом переохлаждения 5Т увеличивается и хтах—х , для 8Т = — 1QK xmax—Хоо = 0.003.
Исследуем влияние скорости охлаждения на распределение состава в твердой фазе. В рассматриваемом режиме выращивания насыщение раствора происходит под влиянием двух факторов: понижения температуры и поступления материала подложки в жидкую фазу (растворение подложки) . При низких скоростях охлаждения основной вклад вносит растворение. Количество поступившего в раствор материала подложки для a = Q.5K/мин примерно в 1.5 раза больше, чем для а = 1.5/Г/мин (см. таб. 2.4). С увеличением а возрастает влияние охлаждения на процесс насыщения раствора. Уменьшается количество поступившего за счет растворения материала и рост начинается при более низкой температуре (см. рис. 2.22). Температура начала роста непосредственно связана с высотой пика концентраций на гетерогранице. Чем выше температура, при которой происходит смена растворения ростом, тем больше неоднородность состава ЭС. 0.218
Определяющее влияние температуры начала роста на величину скачка концентраций на границе подложка - эпитаксиальный слой иллюстрируют результаты расчетов, полученные для различных значений толщины раствора и коэффициентов диффузии в жидкой фазе. На рисунке 2.23 приведено распределение состава твердой фазы для Н1 = 0.15см, 0.2см, 0.25см, скорость охлаждения а = 0.5ііГ/мин, начальное переохлаждение раствора ST = —10К. Естественно, что чем больше толщина жидкой фазы, тем большее количество материала подложки необходимо для ее насыщения (таб. 2.5). Это означает, что с увеличением толщины жидкой фазы растет длительность процесса растворения.
Состав твердой фазы для коэффициентов диффузии D1 = 540 6см2/с, 5 10 5см2/с, 5 10_4см2/с приведен на рисунке 2.24. Толщина раствора Н1 — 0.15см, скорость охлаждения а = 0.5К/мип, начальное переохлаждение 5Т = —10К. Растворение сменяется ростом, когда раствор насыщен вблизи фронта кристаллизации и состав жидкой фазы становится практически однородным (рис. 2.19, 2.20 линия 4). Если диффузия в растворе велика, то перераспределение материала, поступающего из подложки, происходит быстрее и рост начнется при более высокой темпера -24 -20 -16 -12 -8 -4
На рисунке 2.25 показана эволюция состава жидкой фазы на фронте кристаллизации для 5Т = — 2К, — 5К, —10К. Графики представляют собой отрезки прямых, заключенные между изотермами То = 783К и То = 773К. Эти прямые параллельны и, как показали расчеты, их угол наклона практически не зависит от диффузии в жидкой фазе, толщины раствора и скорости охлаждения, а определяется солидусом исходного раствора хо и начальным переохлаждением ST. Как только подложка приводится в контакт с недосыщенным раствором, состав твердой фазы на фронте кристаллизации становится равным хе. При этом хе = хо — АхбТ, где Дх вычисляется по формуле (2.1.3). Если на фазовой плоскости провести прямую, проходящую через точки х0 и Xg, найти точку пересечения этой прямой и изотермы Т = 773К, то получим состав ЭС на поверхности.
Выращивание на подложку менее тугоплавкую, чем солидус исходного раствора. Эволюция состава жидкой фазы на фронте кристаллизации. Распределение CdTe в твердой фазе. Образуется пик. в сторону увеличения или уменьшения содержания CdTe на поверхности твердой фазы (рис. 2.26, 2.27). Если выращивание осуществляется на подложку более тугоплавкую, чем солидус исходного раствора (х, Хо), то содержание CdTe в ЭС будет ниже, чем в подложке (рис. 2.26). Всплеск концентрации CdTe на гетерогранице может возникнуть только в случае слабой диффузии в твердой фазе и при Хоо хо, (хо ф 1) (рис. 2.27).
Граничные условия на фронте кристаллизации
Тестовые расчеты показали, что в случае четырехкомпонентных соединений последовательные алгоритмы, использующие схему 1 или аналогичные ей схемы с другим порядком определения неизвестных, также являются неустойчивыми.
Опираясь на результаты, полученные в [31,67], были предприняты попытки построить для четырехкомпонентных систем численный метод, сочетающий в себе устойчивость совместных алгоритмов и относительную простоту численной реализации последовательных.
Схема 2. о\ = (72 = 1, 7з = О, (J4 = сг5 = сгб = сг7 = 1. В этом случае задача (3.2.2)-(3.2.10) решается следующим образом. На каждом временном слое сначала из уравнений (3.2.2), (3.2.3) с граничными условиями (3.2.5), (3.2.6), (3.2.8), (3.2.9), куда входит c Zn с предыдущего шага по времени, с помощь метода Ньютона определяются c Cd\ с Нд\ Затем найденные значения подставляются в (3.2.7), и полученное таким образом соотношение используется как граничное условие для (3.2.4): c Zn находится из (3.2.4), (3.2.7), (3.2.10).
Схема 3. (7i = аг = сг3 = 1, сг4 = 75 = JQ = 1, 07 = 0. Как и в предыдущем случае, искомые функции на каждом временном слое вычисляются последовательно, сначала совместно определяются &Gd\ с Н9\ а затем досчитывается c Zn . Только в уравнениях для c Cd , с Нд вместо граничного условия (3.2.5) используется (3.2.7), a c Zn определяется из (3.2.4) с граничными условиями (3.2.5), (3.2.10).
Опыт использования схем расщепления типа 2 и 3, в которых сначала совместно определяется распределение любых двух компонент, а затем вычисляется распределение третей, показал, что при любом наборе параметров ап в зависимости от заданных термодинамических параметров и фазовой диаграммы системы не сходятся итерации по методу Ньютона. Таким образом, совместный алгоритм остается единственно возможным.
Схема 4- (?п — 1) п = 1,2,...7. При таком выборе весов разностная схема является полностью неявной и сохраняет нелинейность, существующую в исходной дифференциальной задаче. Для решения системы сеточных уравнений используется итерационный метод Ньютона. Соответствующая система линейных уравнений решается совместно относительно вектора неизвестных, компонентами которого являются все искомые функции. В одномерном случае здесь применяется метод матричной прогонки относительно вектора [c Cd\S-H9\c Zn )T. Как и для трехкомпонентных систем такой подход позволяет построить практически безусловно устойчивый вычислительный алгоритм.
Использование совместных методов предполагает на каждом временном слое решение системы сеточных уравнений с числом неизвестных, равным, в случае трехкомпонентных систем, удвоенному, а в случае четверных соединений, утроенному числу узлов разностной сетки. Свойства матрицы такой системы уравнений, из-за переменных коэффициентов и специфических граничных условий, неизвестны. Поэтому применение итерационных методов, требующих для своей реализации априорной информации об операторах схемы, сопряжено со значительными трудностями. Применение прямых методов в данном случае также малоэффективно, так как требует больших затрат машинного времени. В случае трехкомпонентных соединений [31,67] при решении системы линейных уравнений использовалась специальная процедура переупорядочения неизвестных с последующим применением метода сопряженных градиентов. Попытки обобщить такой подход на случай четверных соединений не дает удовлетворительных результатов из-за плохой сходимости метода сопряженных градиентов.
