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

Введение.

1 Методы численного решения уравнений динамики вязкой несжимаемой жидкости.

1.1 Схемы с искусственной дисперсией для уравнения конвективной диффузии.

1.1.1 Линейная модельная задача.

1.1.2 Анализ устойчивости схем с искусственной дисперсией.

1.1.3 Разностная схема с искусственной дисперсией для одномерного уравнения конвективной диффузии.

1.1.4 Разностная схема с искусственной дисперсией для двумерного уравнения конвективной диффузии.

1.1.5 Сравнение разностных схем.

1.2 Матричный алгоритм решения уравнений Навье - Сток

1.2.1 Общие замечания.

1.2.2 Исследование устойчивости модельной задачи.

1.2.3 Линейный матричный алгоритм.

1.3 Основные результаты, полученные в главе 1.

1.4 Рисунки к главе 1.

2 Исследование устойчивости алгоритмов численного решения систем параболических уравнений с нетрадиционными граничными условиями.

2.1 Жидкофазовая эпитаксия трёхкомпонентных соединений. Диффузионная модель.

2.1.1 Линеаризация задачи.

2.1.2 Анализ одномерной дифференциальной модели.

2.1.3 Устойчивость разностных схем для модельной задачи.

2.2 Задача о фазовом переходе в бинарной системе. Одномерная модель.

2.2.1 Формулировка модельной задачи.

2.2.2 Исследование устойчивости разностных схем для модельной задачи о кристаллизации бинарных соединений.

2.3 Основные результаты, полученные в главе 2.

Математическое моделирование нестационарных процессов кристаллизации бинарного раствора.

3.1 Описание технологического процесса.

3.2 Математическая постановка задачи.

3.3 Метод численного решения.

3.3.1 Замена переменных.

3.3.2 Дифференциальные операторы в новых переменных.

3.4 Разностная схема.

3.5 Описание алгоритмов численного решения разностной задачи.

3.6 Сравнение методов.

3.7 Моделирование процесса кристаллизации бинарного соединения в ампуле.

3.8 Основные результаты, полученные в главе 3.

3.9 Рисунки к главе 3.

Численное исследование процесса получения двойных полупроводниковых соединений методом жидкофазо-вой эпитаксии.

4.1 Математическая постановка задачи.

4.2 Выращивание эпитаксиальных слоев фосфида галлия в вертикальных системах.

4.3 Определение скорости травления подложек фосфида галлия в ненасыщенном расплаве Ga — GaP.

4.4 Выращивание эпитаксиальных слоёв в "градиенте температуры".

4.5 Выращивание эпитаксиальных слоёв в горизонтальных системах.

4.6 Основные результаты, полученные в главе 4.

4.7 Рисунки к главе 4.

5 Математическое моделирование эпитаксиального выращивания твёрдых растворов Сс1уН^1-уТе из жидкой фазы.

5.1 Математическая постановка задачи.

5.2 Численное исследование процесса эпитаксиального выращивания СйуНд\уТе из растворов - расплавов на основе Те.

5.2.1 Диффузионное приближение, область применимости.

5.2.2 Влияние структуры конвективного движения на однородность распределения состава и толщину эпитаксиальных слоев.

5.3 Выбор режимов получения слоёв СйуНд\уТе с гладкой поверхностью и однородным распределением состава.

5.4 Основные результаты, полученные в главе 5.

5.5 Рисунки к главе 5.

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

Новый импульс развитию вычислительной динамики вязкой жидкости дали проблемы, возникающие в полупроводниковой технологии, в частности задачи, связанные с получением монокристаллов и полупроводниковых структур из жидкой фазы, математическое моделирование которых интенсивно проводится в последние десятилетия. Высокие требования, предъявляемые современным приборостроением к свойствам и качеству полупроводниковых структур и материалов, могут быть удовлетворены лишь на основе детального изучения физических процессов, происходящих при их получении, и тщательной отработки самого технологического процесса. Однако натурные эксперименты в этой области сложны и трудоёмки, требуют больших энергетических затрат, в них используются дорогостоящие, зачастую ядовитые материалы, а длительность одного эксперимента может достигать нескольких недель. Поэтому использование методов математического моделирования в исследованиях по разработке, усовершенствованию и оптимизации технологических режимов и установок, для анализа результатов физических экспериментов является необходимым этапом на пути создания способов получения полупроводниковых материалов с заданными свойствами.

Диссертация посвящена разработке и исследованию методов численного решения нестационарных задач конвективного тепло - мас-сопереноса в системах с фазовыми переходами. Алгоритмы ориентированы на изучение процессов кристаллизации многокомпонентных (двойных и тройных) соединений и применяются для анализа двух технологических процессов: метода жидкофазовой эпитаксии и метода вертикальной направленной кристаллизации.

Метод жидкофазовой эпитаксии (ЖФЭ) - это способ выращивания монокристаллических слоёв из раствора - расплава на подложке определённой кристаллографической ориентации [1]. В простейшем случае получения плёнок двойного полупроводникового соединения АВ подложка из этого материала приводится в контакт с раствором компонента В в расплаве А. Если при температуре контакта (Го) раствор - расплав является насыщенным, то жидкая и твёрдая фазы находятся в квазиравновесии. Для того чтобы вызвать процесс кристаллизации, температура системы расплав - подложка понижается с заданной скоростью. С понижением температуры падает растворимость компонента В в А, раствор становится пересыщенным и из него на подложке кристаллизуется соединение АВ. В процессе роста эпитаксиального слоя на границе раздела фаз сохраняется квазиравновесие между раствором - расплавом и эпи-таксиальным слоем, вместе с тем основной объём расплава остаётся пересыщенным. Квазиравновесный режим кристаллизации означает, что состав жидкой фазы на границе расплав - кристалл изменяется по кривой растворимости. В случае получения тройных полупроводниковых соединений состав жидкой и твёрдой фаз определяется фазовой диаграммой рассматриваемой системы. Перенос растворённых компонентов к фронту кристаллизации осуществляется механизмом конвекции и диффузии.

Гидродинамические процессы, протекающие в растворе - расплаве, оказывают существенное влияние на свойства эпитаксиаль-ных слоёв и их качество, определяя характер поступления растворённых компонентов к поверхности роста и тем самым скорость роста, равномерность её распределения вдоль поверхности эпитаксиального слоя (ЭС), а в многокомпонентных системах - состав кристаллизующегося материала. Конкретные особенности массопе-реноса в растворе - расплаве зависят от условий технологического процесса: объёма расплава, скорости его охлаждения, положения подложек (вертикальное, горизонтальное, над расплавом или под ним), величины и направления температурных градиентов и т.д.

Метод вертикальной направленной кристаллизации [2] используется для получения объёмных монокристаллов из жидкой фазы.

Цилиндрическая ампула, на дне которой находится затравка из кристаллизуемого материала, заполняется расплавом (раствором) соответствующего соединения. Ампула помещается в печь, где специальным образом создаётся тепловое поле. Внутри печи ампула перемещается из области высоких температур в область более низких и постепенно охлаждается. В результате расплав, находящийся в ампуле, кристаллизуется.

При получении чистых веществ или материалов, легированных небольшим количеством примесей, граница расплав - кристалл проходит по изотерме со значением температуры, равным температуре фазового перехода. В данной работе рассматриваются процессы кристаллизации двойных соединений, температура фазового перехода которых зависит от состава жидкой фазы и определяется фазовой диаграммой системы. Условия кристаллизации, скорость роста, граница фазового перехода, состав твёрдой фазы формируются в результате сложного взаимодействия процессов теплопереноса в кристалле и ампуле, теплообмена между ампулой и окружающей средой, естественной конвекции в расплаве, вызванной неоднород-ностями в распределении температуры и концентрации и т.д.

В диссертации численное исследование метода жидкофазовой эпи-таксии и метода вертикальной направленной кристаллизации проводится на основе решения двумерных нестационарных уравнений тепловой и концентрационной конвекции в приближении Буссине-ска [3].

Общая характеристика рассматриваемого класса задач. Эффективность применения методов математического моделирования для изучения явлений тепло - массопереноса в конкретных системах и установках в значительной степени определяется качеством используемых алгоритмов. В полной мере это относится к задачам полупроводниковой технологии, в которых необходимо проводить многопараметрические исследования, связанные с анализом результатов большого числа расчётов, сопоставлением их с данными натурных экспериментов. С точки зрения практического применения результатов математического моделирования, методы должны быть ориентированы на использование сравнительно маломощных вычислительных средств, непосредственно имеющихся в распоряжении специалистов - технологов, и обеспечивать получение надёжных результатов при разумных затратах машинного времени на расчёт каждого варианта. Для того чтобы удовлетворить этим требованиям, вычислительный алгоритм необходимо конструировать с учётом конкретных особенностей изучаемых процессов и явлений. Важное значение имеет определение реальной области применимости разработанных численных методов, того диапазона изменения физических параметров задачи, в котором алгоритмы гарантируют точность получаемых результатов.

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

V 1 о лов с числом Прандтля Рг — — и 10 -Ь 10 и числом Шмидта Л

5с = — « 5 х 101 -г-102. (V, А, и Б - соответственно коэффициенты кинематической вязкости, температуропроводности расплава и диффузии растворённых компонентов.) Это означает, что масштабы времени, связанные с диссипативными процессами: теплопроь2 ь2 ь2 /т водностью £д = —, вязкостью — — и диффузиеи ¿о = — (Ь ли и характерный линейный масштаб задачи) отличаются между собой на один - два порядка: Рг = —, 5с = —, < < Ьр.

Ещё один масштаб времени (иес) в этих задачах определяется технологическими режимами, которые должны, например, обеспечить необходимую толщину эпитаксиального слоя или максимальный объём монокристалла с заданными свойствами. Длительность процесса ¿¿ес непосредственно связана со скоростью выращивания кристаллов Vcr : ttec = ——, где l - толщина ЭС или длина кристалла. При моделировании нестационарных процессов ttec задаёт отрезок времени, на котором необходимо находить решение задачи. Обычно он на порядок и более превосходит время вязкой диссипации и в несколько раз больше диффузионных времён, а ~ 106. t\

Трудности численного исследования нестационарных процессов кристаллизации в первую очередь связаны с широким диапазоном изменения временных масштабов, присутствующих в задаче. Используемые алгоритмы должны обеспечивать необходимую точность расчёта быстрых процессов и, вместе с тем, не требовать слишком больших затрат машинного времени на расчёт всей задачи в целом. Важно, чтобы выбор шага по времени определялся спецификой физических процессов, протекающих в системе, а не диктовался внутренними свойствами вычислительной процедуры. В частности это означает, что для численного исседования методов получения полупроводниковых материалов из жидкой фазы необходимо использовать такие разностные алгоритмы, в которых величина шага по времени жёстко не связана с пространственным масштабом разностной сетки. Прежде всего это относится к алгоритмам решения уравнений движения.

Традиционные методы численного решения уравнений Навье -Стокса в переменных "функция тока (ф) - вихрь (и>) " [4], основанные на раздельном на каждом временном слое определении ш и ф, обладают ограничением на величину максимально допустимого в расчётах шага по времени: т < г* « /гтогп = min uL пТ7 nv - шаги разностной сетки по пространству, Re = — - число v

Рейнольдса.

Требование г < т* является оправданным при моделировании интенсивной вынужденной конвекции (Re 1), когда скорости движения жидкости велики и решение нужно определить на сравнительно небольшом промежутке времени t ~ tv. При Re 1 размер минимального шага по пространству выбирается исходя из с Ь ширины динамического пограничного слоя д ~ , так что шаг

V В,е по времени оказывается не слишком мелким и хорошо согласуется с внутренней структурой течения. Однако и в расчётах медленных ("ползущих") течений с малыми значениями числа Рейнольдса или в случае естественной конвекции с небольшими значениями диффузионного и теплового числа Грасгофа, когда поле скоростей изменяется медленно и расчёт уравнений движения можно проводить на грубых временных сетках, величина максимально допустимого в расчётах шага по времени по - прежнему ограничена условием г « /¿тот- Оно оно уже не связано с характером гидродинамических процессов и обусловлено исключительно свойствами численного метода.

Ограничение т < т* становится особенно обременительным при моделировании концентрационной конвекции в средах с большими значениями числа Шмидта, когда определяющим размером для величины шага по пространству является ширина диффузионно

А 1 ^ го пограничного слоя (— « ) и отрезок времени, на котором д \/5с необходимо определить решение, велик: £ > = 5с Именно такая ситуация возникает в нестационарных задачах кристаллизации. Ограничение т ~ требует неоправданно большого числа шагов по времени и, следовательно, больших затрат машинного времени для определения искомых функций на всём отрезке Ь < представляющем практический интерес. Интенсивность конвективного движения в расплаве невелика: реальные значения теплового д[ЗтЬ3АТ д(ЗсЬ'АС иг? =--- и диффузионного игр =--- числа 1 расго

V1 и1 фа лежат в диапазоне 104-г 107. Поэтому с физической точки зрения уравнения движения можно решать крупным шагом по времени и тем самым сократить общее время расчёта задачи. ( ДТ, АС

- характерные перепады температуры и концентрации, (Зт, (Зс

- коэффициенты теплового и концентрационного расширения, д -ускорение свободного падения.)

Важную роль в конструировании алгоритмов численного исследования задач кристаллизации играет выбор метода решения уравнения конвективной диффузии. Предполагая, что характерная скорость движения жидкости (по) имеет порядок 0(л/дЬ[ЗАТ) [5, стр.63, щЬ

341], диффузионное число Пекле Реп = -ц- можно оценить величиной л/ОгЗс >> 1, т.е. массоперенос в расплаве описывается уравнением конвективной диффузии с преобладающим влиянием конвекции.

Хорошо известны трудности, возникающие при численном решении таких уравнений [7] - [10, т.1] : замена конвективных членов направленными разностями первого порядка вносит в ошибку аппроксимации слагаемые, которые по величине могут превосходить члены, описывающие реальный процесс диффузии; симметричная запись конвективных членов с помощью центральных разностей приводит в областях резкого изменения решений к возникновению "нефизических колебаний", обусловленных внутренними свойствами схемы.

Необходимую точность центральные и направленные разности гарантируют в случае, когда локальное значение сеточного числа ик

Пекле Рв}1 = — < 2, т.е. на сетках с шагом по пространству г « — и числом узлов по каждому направлению N ~ \fGrSc. \fGrSc

Здесь, как и в случае уравнений Навье - Стокса, требования к параметрам разностной сетки связаны не со свойствами изучаемого физического процесса, а с особенностями численного метода. Величина шага по пространству в рассматриваемых задачах определяется шириной диффузионного пограничного слоя 6с ~ 4/ fGr\JБс поэтому надёжные результаты могут быть получены на сетках с числом узлов N ~ \fGr\TSc <С N. Применение алгоритмов, обеспечивающих высокую точность расчётов при физически оправданных ограничениях на параметры сетки, позволяет на порядок и более сократить число узлов сетки по каждому направлению. Точность на грубых сетках в первую очередь необходима в расчётах задач кристаллизации многокомпонентных растворов, когда поле концентрации влияет на основные характеристики процесса: положение границы фазового перехода и состав твёрдой фазы.

Ещё одной особенностью рассматриваемого класса задач является наличие граничных условий специального вида. Они возникают при математическом описании равновесного процесса кристаллизации многокомпонентных растворов (расплавов) [11, т.1, стр.388]. На границе раздела фаз температура и концентрация должны удовлетворять не только балансным соотношениям, вытекающим из законов сохранения энергии и массы, но и уравнению ликвидуса фазовой диаграммы системы: довольно сложному, как правило, нелинейному соотношению, связывающему равновесные концентрации всех компонентов и температуру.

Для моделирования процессов кристаллизации многокомпонентных соединений обычно используются алгоритмы, в которых аппроксимация граничных условий позволяет на каждом временном слое решать уравнения переноса тепла и/или массы для каждого из компонентов последовательно, что достигается за счёт введения в аппроксимацию граничных условий элемента явности. На границе значение первой из искомых функций на новом временном слое выражается через значения с предыдущего слоя остальных функций, входящих в соответствующее граничное условие. Затем оставшиеся уравнения решаются последовательно с использованием в граничных условиях значений функций, уже найденных на текущем слое (процедура типа Зейделя). Эти методы экономичны и просты в реализации, однако, обеспечивая получение правильных результатов в одном диапазоне параметров (начальные данные, фазовая диаграмма, термодинамические константы и т.д.), в другом - они оказываются неустойчивыми или дают правдоподобные, но неверные результаты [12],[13]. Проблема состоит в том, как зная фазовую диаграмму и термодинамические свойства системы, выбрать способ численной реализации граничных условий, обеспечивающий устойчивость вычислительной процедуры в рассматриваемом диапазоне изменения параметров.

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

• разностных схем, гарантирующих на сетках со сравнительно небольшим числом узлов по пространству высокую точность решения уравнений конвективной диффузии с малым параметром при старшей производной;

• метода решения уравнений Навье - Стокса, позволяющего использовать грубые сетки по времени;

• надёжного способа численной реализации условий на границе раздела фаз.

Рассмотрим некоторые из применяемых на практике способов конструирования таких алгоритмов.

О методах решения уравнений типа конвективной диффузии. Аппроксимация конвективных членов с первым порядком точности порождает в ошибке аппроксимации слагаемые со вторыми пространственными производными, которые на грубых сетках вместо реальных диффузионных членов определяют диссипацию в дискретной среде. Поэтому обычно разностные алгоритмы для решения уравнений с преобладающим влиянием конвекции строятся в классе схем более высокого порядка точности. Их можно интерпретировать как схемы с центрально - разностной аппроксимацией конвективных членов, в которые введены дополнительные слагаемые - сеточные аналоги пространственных производных третьего, четвёртого порядка или смешанные производные. Они играют роль регуляри-заторов, позволяющих влиять на дисперсионные и диссипативные свойства алгоритмов и их устойчивость.

Общие принципы регуляризации разностных схем разработаны А.А.Самарским (см., например, [7, стр.345], [11, т.1, стр.280]). В конкретных случаях построение регуляризатора проводится на основе различных идей и подходов. Отбор и оценка полученных алгоритмов также осуществляются с помощью различных критериев: выполнение принципа максимума для соответствующего сеточного оператора [14], [15], порядок аппроксимации схемы и её шаблон

7, стр.147], [16] - [19], поведение корней характеристического уравнения [20] - [22] или поведение Фурье - компонент приближённого решения [23], [10, т.1 стр.373, 406], свойств матрицы разностного оператора и т.д.

Широко используются алгоритмы, в которых регуляризатор конструируется исходя из непосредственного анализа локальных свойств решения [24] - [32] и внутри расчётной области в зависимости от поведения решения происходит переключение с одной схемы на другую.

В алгоритмах второго порядка точности подавление конвективных членов за счёт увеличения схемной диффузии в областях с большими значениями локального числа Пекле может осуществляться автоматически. Из большого количества работ в этом направлении отметим лишь недавнюю статью [33] и примыкающую к ней по смыслу работу [34], где предложены схемы типа экспоненциальных. Асимптотически при больших значениях сеточного числа Пекле Ре^ они переходят в схемы с аппроксимацией конвективных членов направленными разностями, а при стремлении Ред/г к нулю конвективные члены в разностном уравнении переноса пропадают.

Регуляризация схем с центральными разностями с помощью третьих или четвёртых пространственных производных проводилась в работах [35] - [38], где дополнительные слагаемые, имеющие порядок 0(Н2) или 0(/г3) выбирались таким образом, чтобы на неравномерной сетке и при переменной скорости регуляризаторы были знакоопределёнными операторами и обеспечивали безусловную устойчивость схем в норме ¿2- Обычные требования к свойствам матрицы сеточного оператора, поведению частных решений и т.д. также сохранялись.

Для изучения свойств сеточных алгоритмов, построения регуля-ризаторов, их классификации широко используется метод дифференциального приближения [39]. В зарубежной литературе аналогичный способ исследования разностных схем получил распространение под названием метода модифицированного уравнения. (См., например, [41], [10, т.1].) Анализ структуры ошибки аппроксимации, при котором производные чётного и нечётного порядков связываются соответственно с диссипативными и дисперсионными характеристиками разностного алгоритма, даёт наглядную физическую интерпретацию поведения решения дискретной задачи. Введение в разностную схему диссипативных или антидисперсионных добавок, компенсирующих в той или иной форме влияние старших членов в ошибке аппроксимации, позволяет приблизить свойства разностных решений к дифференциальному случаю, при этом формальный порядок аппроксимации метода может и не измениться.

Систематическое применение метода дифференциального приближения связано в первую очередь с задачами газовой динамики [39], в частности с его помощью были построены и изучены схемы с искусственной дисперсией для нелинейного уравнения переноса, одномерных уравнений газовой и магнитной газовой динамики [42]

44]. Для конструирования алгоритмов решения уравнений конвективной диффузии с малым параметром при. старшей производной метод дифференциального приближения использовался в работах

45] - [48].

На основе анализа модифицированного уравнения разностные схемы для линейного уравнения переноса и уравнения конвективной диффузии с постоянной скоростью изучались также в [49] - [53]. В этих работах конструирование схем осуществлялось в рамках формального подхода. Аппроксимация уравнений записывалась в виде суммы с неопределёнными коэффициентами разностных отношений по времени и пространству (направленных и центральных разностей). Коэффициенты затем определялись таким образом, чтобы в рассматриваемом классе схем в том или ином смысле оптимизировать диссипативные и дисперсионные свойства алгоритма. В результате были получены многие из известных ранее схем [50]. Новые же методы, например третьего порядка, конструировались главным образом без расширения шаблона, как линейная комбинация схем второго порядка.

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

Анализ устойчивости полученных схем проводится, как правило, с помощью метода Неймана [10, т.1],[41].

В данной работе для численного решения уравнений массопере-носа используется разностная схема с искусственной дисперсией [46] - [48], в которой к центрально - разностной аппроксимации конвективных членов добавлено слагаемое порядка О (/г2) с направленными третьими разностями, компенсирующее дисперсионное влияние старшего члена в ошибке аппроксимации. При постоянной скорости в одномерном приближении у схемы отсутствует схемная дисперсия третьего порядка ( в случае, когда концентрация относится к узлам сетки и ошибка аппроксимации определяется в тех же узлах). В двумерных задачах с переменной скоростью введённый в схему регу-ляризатор с третьими пространственными производными по каждому направлению не уничтожает полностью член порядка О (к2) в ошибке аппроксимации, однако он позволяет существенно улучшить дисперсионные свойства алгоритма.

Методические расчёты показывают, что схемы с искусственной дисперсией обеспечивают получение надёжных результатов на сетках с шагом по пространству /г « у== [46], [47]. На таких же сравнительно грубых сетках схемы с искусственной дисперсией дают возможность моделировать нелинейные процессы массопереноса, когда скорость в уравнении конвективной диффузии зависит от концентрации, как это имеет место при электрофоретическом

Ьи

разделении биосмесей [55], [56]. Значение числа Пекле Ре = в задачах электрофореза изменяется в диапазоне 104 -Ь 106 [57], [58].

О методах решения уравнений Навъе - Стокса. Рассмотрим теперь некоторые из алгоритмов численного решения уравнений движения вязкой несжимаемой жидкости. В данной работе моделирование процессов конвективного тепло - массопереноса проводится на основе уравнений Навье - Стокса, записанных в переременных "вихрь и), функция тока 1р." Такая форма записи широко используется при решении задач в двумерном приближении [4].

Работы по созданию разностных схем для уравнений параболического типа с малым параметром при старшей производной в значительной степени стимулировались интересом к проблемам, связанным с моделированием движения жидкости при больших значениях числа Рейнольдса. Как и в случае уравнений конвективной диффузии с преобладающим влиянием конвекции, усилия были направлены на то, чтобы построить схемы порядка выше первого, у которых отсутствовали бы "нефизические" колебания в областях с резким изменением решения, а матрица разностного оператора обладала бы свойством диагонального преобладания [4], [10, т.1], [15], [59] - [61].

В частности, компактные разностные схемы порядка 0(/г4) и выше были построены в [18], [19], [62]. Используя стандартную процедуру повышения порядка точности [7, с.249], схемы конструировались на девятиточечном шаблоне. (Уравнение в узле (г, ]) содержит значения искомой функции в этой точке и в восьми её ближайших соседях.) Выполнение свойства диагонального преобладания достигается за счёт использования направленных разностей в членах, введённых в схему для повышения порядка аппроксимации.

Важным этапом в развитии методов решения уравнений Навье

- Стокса явилась работа [63]. В дифференциальном случае члены дш дф ди; дф переноса V = ——-----—— не дают вклада в баланс кинетической ох оу оу ох

1 и? энергии £ = -(и2 + V2) и в баланс величины Б — — (энстропии). В и ¿А работе [63] было предложено семейство схем, сохраняющих на сетке величины £ и S. Схемы, в которых аппроксимация членов переноса не даёт вклада в баланс кинетической энергии, позднее получили название энергетически нейтральных [65].

Как показано в [66], нарушение энергетической нейтральности делает невозможным моделирование гидродинамической неустойчивости: решение быстро возрастает и теряет физический смысл. При изучении других процессов энергетические свойства разностных схем могут не играть столь решающей роли, однако аппроксимация А.Аракавы [63] позволяет сохранить в дискретной модели балансные соотношения, присущие исходной дифференциальной задаче, что важно для построения надёжной вычислительной процедуры [9].

Дальнейшие исследования [67] позволили выделить среди предложенных в [63] схем, основанных, как известно, на аппроксимации конвективных членов центральными разностями, энергетически нейтральную разностную схему, позволяющую на грубой сетке с хорошей точностью моделировать движение жидкости при достаточно больших значениях числа Рейнольдса (Яе ~ 103). В данной работе конвективные члены аппроксимируются с помощью этой схемы.

Формулировка краевой задачи для уравнений Навье - Стокса в переменных "вихрь, функция тока" формально не содержит граничных условий для ш. Физически согласованные условия можно вывести для ф и —^—. Вытекающее из них условие на со носит нелокаль-оп

I ^ ныи характер: пусть ф\т = дп 0 и £ - функция удовлетворяюг щая уравнению Лапласа в области S, тогда справедливо равенство f u^dS — 0, т.е. вихрь ортогонален в смысле L2 подпространству s гармонических функций. Интегральное условие не сводится к локальному заданию вихря на границе [68]. По существу попыткам разрешить это противоречие в рамках методов, основанных на последовательном на каждом временном слое решении уравнений для (V и ф, и посвящены работы, в которых строятся и изучаются различные способы вычисления вихря на твёрдой стенке [4], [8, с.214], [10, т.2, с.442], [69] - [73].

Способ получения граничних условий обычно состоит в разложении функции тока в ряд вблизи границы. В результате с той или иной точностью определяется связь между функцией тока и вихрем: (¿|г = f(4>). При раздельном решении значение ¿D|r на верхнем временном слое вычисляется по значениям ф с предыдущего слоя по времени или с предыдущей итерации. Уравнение Пуассона для определения функции тока решается с граничным условием ф|г =0. В ряде работ [74], [4, с.190] выполнение на твёрдой стенке условия 0 достигается с помощью непосредственного г ^ прилипания -т— on подправления функции тока в окрестности границы. При этом сами уравнения в некоторой приграничной полосе могут не выполняться. Для устранения невязки в нестационарных задачах необходимо использовать на каждом временном слое внутренний итерационный процесс.

Методы, основанные на последовательном решении уравнений для функции тока и вихря, обычно обладают ограничением на величину максимально допустимого в расчётах шага по времени т < 3 г* = 7 min (/4, hy)Re, 7 < - [75] - [81]. (7 определяется конкретным 2 видом граничных условий, их точностью, согласованностью с разностной схемой, используемой во внутренних узлах и т.д. [8, с.214], [72], [73].) Если г > г*, то неустойчивость первоначально проявляется в поведении вихря. Она зарождается в приграничных узлах, затем распространяется внутрь области [75], [78].

Одной из основных причин, ограничивающих величину шага по времени в алгоритмах с последовательным решением уравнений переноса вихря и Пуассона, является локальное задание граничных условий для икоторое не обеспечивает выполнение интегральных соотношений, справедливых в дифференциальном случае. Вычисление вихря на границе по значениям функции тока с предыдущего слоя фактически означает, что, вместо разностного аналога соотношения J toÇdS = 0, для сеточных функций выполнено равенство — ri*1, где ~ скалярное произведение в пространстве сеточных функций, ^ - разностная гармоническая функция, г -шаг сетки по времени, F - невязка.

Если î[j\y = 0, w|r определяется по формуле Тома [82], что, как показано в [65], необходимо для получения консервативного алгоритма, и для наглядности £/,, = 1, то на равномерной прямоугольной h сетке с точностью до членов порядка 0(т) F = ]Г) ~r~('ll;t{x2,yj) + j=2 hх

M~l hx t{xM~uVj)) + I] ir-(i>t(xi,y2) + 1)) (hx,hy - шаги разi= 2 "y ностной сетки, (а^, у у) - её узлы, г = 1,., М, = 1, . . , И).

Попытки повысить устойчивость вычислительной процедуры, вынося в граничных условиях для вихря функцию ф на верхний временной слой и оставаясь при этом в рамках алгоритмов с раздельным решением уравнений для ш и ф, предпринимались многими авторами. (См., например, [79], [83] - [85].)

В [79] для определения функции тока на верхнем слое решается уравнение переноса вихря, в котором в производной по времени со на верхнем слое заменяется на — Дф, с верхнего же слоя берётся функция тока в конвективных членах и граничных условиях Тома. Входящие в это уравнение значения вихря относятся к нижнему слою. Расчёт вихря производится из того же уравнения, записанного в обычной форме с использованием уже найденных значений функции тока.

В [83] ослабление влияния невязки ^ на устойчивость численного метода достигается за счёт введения в уравнения для определения функции тока дополнительного слагаемого, отличного от нуля вблизи границы. Оно пропорционально ^ и имеет порядок О (г2).

В [68] условие ортогональности вихря подпространству сеточных гармонических функций непосредственно участвует в определении со.

В алгоритмах, основанных на совместном решении уравнения переноса вихря и уравнения Пуассона, когда на каждом временном слое система уравнений Навье - Стокса рассматривается как одно матричное уравнение относительно двухкомпонентного вектора м необходимости в специальном граничном условии для ш не

Ф/ дф возникает: как и в дифференциальном случае условии на ф и —— оп достаточно для получения решения. Матричное уравнение на гра п дФ нице непосредственно выписывается с учетом того, что ф г и —— оп г тт ЭФ заданы. При согласованном выборе аппроксимации производной ——

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

Одна из первых работ, где сеточные уравнения Навье - Сток-са рассматривались как связанная система уравнений относительно вектора ( ), принадлежит Н.И.Булееву и Г.И.Тимухину [86].

Ф)

Однако использовавшийся ими итерационный метод решения системы линейных матричных уравнений [87] не обладал достаточной скоростью сходимости, что не позволило получить заметный выигрыш в величине шага по времени по сравнению с раздельными алгоритмами.

Автором данной работы был предложен метод совместного решения нестационарных уравнений Навье - Стокса, оказавшийся практически безусловно устойчивым [75] - [78]. С его помощью моделирование нестационарных процессов естественной конвекции с умеренными значениями числа Рэлея (Да) можно проводить шагом по времени т « 102 -Ь 103 тш (/г2, Ь?х). Хотя число арифметических действий на временном слое в этом методе существенно больше, чем в последовательных алгоритмах, благодаря использованию крупных шагов по времени удаётся существенно сократить затраты машинного времени на расчёт задачи в целом.

Алгоритм основан на неявной по ф иш аппроксимации уравнения переноса вихря. Граничные условия для со также содержат ф с верхнего слоя по времени. Нелинейная система уравнений решается методом Ньютона. Возникающая на ньютоновских итерациях система линейных уравнений относительно приращений 5фц рассматривается как совокупность матричных уравнений относительно векторов 5Хц = ( %3 ), так что в каждой точке расчёт

КНгз) ной области задаётся уравнение, коэффициенты которого являются матрицами второго порядка. Вектор приращений находится с помощью итерационного процесса [88], [89].

Алгоритмы, в которых уравнения Навье - Стокса решаются совал местно относительно вектора Хц = будем называть матрич

Фг3/ ными или, следуя зарубежной терминологии, связанными, совместными (coupled algorithms). [59], [10, т.II].

Связанный алгоритм для решения стационарных уравнений Навье - Стокса был построен в [59]. Там вектор неизвестных определялся с помощью одного из вариантов метода неполной факторизации [90], [10, т.1, стр.258], [91]. Дальнейшие исследования авторов показали [60], что алгоритм [59] мало пригоден для моделирования нестационарных задач: приемлемая точность в методе [59] может быть получена лишь при использовании шага по времени т <С /г2Де. Причиной возникновения столь жёсткого ограничения на величину шага по времени при совместном решении уравнений для со и гр является способ решения уравнений относительно векда неполной факторизации позволила в расчётах модельных задач ослабить ограничение на т, однако как поведёт себя этот вариант метода в реальных задачах остаётся неясным.

Матричные алгоритмы вообще очень чувствительны к выбору способа решения сеточных уравнений. В последовательных алгоритмах со на верхнем слое может определяться с помощью различных экономичных итерационных процедур [92] - [96]. Если схема для уравнения переноса вихря обладает свойством диагонального преобладания, то, как правило, итерации хорошо сходятся и не возникает дополнительных, помимо г < г* « ограничений на шаг по времени [80], [81]. Решение уравнения Пуассона также не вызывает трудностей [4].

В матричных алгоритмах, когда уравнения Навье - Стокса решаются как связанная система, сходимость итерационных методов резко ухудшается. Не раз отмечалось [75], [60], [79], [97], что итерационные процессы типа Зейделя, переменных направлений, локальной релаксации и т.д. в случае, когда их проводят относительно вектора удовлетворительно сходятся лишь на сетках с небольшим числом узлов. В наших опытах по использованию итерационного метода ОКГОМШ [95], [96] для совместного решения уравнений для со и яр искомые функции удалось получить лишь

Предложенная в [60] модификация метона сетке 11x11, чисто неявная процедура [90] в [59] применялась на сетке 17x17; число итераций изменялось в диапазоне 50 -Ь 100. В нестационарных задачах ограничения, которые накладываются на величину шага по времени условиями сходимости итерационной процедуры, могут оказаться ещё более жёсткими, чем те, что имеют место в последовательных алгоритмах [75], [60].

При совместном решении уравнений для ш и ф хорошо зарекомендовал себя метод [89] - матричный аналог предложенного в [88] "а; — /3" итерационного процесса. В сочетании с нелинейным матричным алгоритмом этот метод использовался при моделировании жидкофазовой эпитаксии двойных полупроводниковых соединений [98] - [102]. Расчёты проводилось на сетках с числом узлов по пространству М х N « 20 х 100 шагом по времени т « (102 -г- 103)/г^п, который выбирался не из соображений устойчивости вычислительной процедуры, а из соображений точности получаемых результатов. Метод решения сеточных уравнений не накладывал каких - либо специальных ограничений на величину т. Однако, по сравнению, скажем, с градиентными методами, "си — /3" алгоритм требует для своей реализации больших затрат оперативной памяти и большего числа арифметических действий на каждой итерации. Специальных исследований требует также возможность обобщения этого алгоритма на случай схем, записанных на расширенном шаблоне.

И.В. Фрязиновым был построен метод совместного решения уравнений для ш и ф, в котором матрица сеточного оператора, действующего на вектор является суммой симметричной положительно определённой матрицы и кососимметричной части [103]. Такая структура матрицы, полученная за счёт введения в уравнение Пуассона специальным образом подобранного регуляризатора, позволяет доказать разрешимость сеточной задачи на каждом временном слое. Положительная определённость оператора обеспечивает сходимость итерационного процесса типа ОКГОМШ [95], [96]. Но она оказалась медленной и в [103], [136] для решения системы сеточных уравнений использовался прямой метод с переупорядочиванием

104], [105].

В данной работе для решения уравнений Навье - Стокса используются матричные алгоритмы [76] и [106], последний в дальнейшем будем называть линейным. В отличие от нелинейного алгоритма [76], в методе [106] функция тока в конвективных членах, записанных по формулам А.Аракавы, берётся с предыдущего шага по времени, поэтому система сеточных уравнений на слое становится линейной. Регуляризация конвективных членов осуществляется с помощью введения в схему дополнительных слагаемых порядка О (/г,3) с четвёртыми разностями [35], [37]; для совместного решения уравнений движения использовался прямой метод [104] в сочетании с процедурой переупорядочивания [105].

Линейный метод совместного решения уравнений Навье - Стокса [106], так же как и матричный алгоритм [75] - [78], оказался практически безусловно устойчивым. Прямой метод позволяет избежать каких - либо ограничений на величину шага по времени, связанных со сходимостью итерационных процессов при решении сеточных уравнений. Основной же недостаток прямых методов: резкое увеличение затрат оперативной памяти и машинного времени с ростом числа узлов разностной сетки, удаётся компенсировать благодаря применению более точных разностных схем и учёту специфики рассматриваемых задач, в которых движение жидкости изучается в областях с отношением сторон 1 : 5 -Ь 1 : 20. Это позволяет использовать сетки со сравнительно небольшим числом узлов по одному из направлений: т = min (М, N) < 25. В таких случаях линейный метод в 2 - 5 раз уступает по затратам машинного времени на слое раздельным алгоритмам, а выигрыш в величине шага по времени может составлять 100 и более раз.

О численной реализации условий на границе раздела фаз. Рассмотрим сначала процесс получения тройных полупроводниковых соединений методом жидкофазовой эпитаксии. В рамках предположения о равновесном характере кристаллизации состав жидкой и твёрдой фаз на границе раздела удовлетворяет диаграмме состояния системы. Кроме того, для каждого компонента на границе расплав - кристалл должен выполняться закон сохранения массы. Он связывает поток растворённого компонента со скоростью движения фронта кристаллизации. (В этих задачах диффузионный перенос в эпитаксиальном слое не учитывается; зависимость температуры от времени считается заданной. [1], [107].) Таким образом массоперенос в многокомпонентной среде описывается системой уравнений с граничными условиями специального вида и распределение каждого из компонентов в растворе - расплаве не может определяться независимо: на границе жидкое - твёрдое связаны значения концентраций всех растворённых компонентов и их потоки. В дифференциальном случае скорость движения фронта кристаллизации может определяться по потоку любого из растворённых компонентов.

Аналогичная ситуация возникает и при моделировании роста объёмных монокристаллов из бинарного расплава. Здесь температура фазового перехода зависит от состава жидкой фазы, а поток растворённого компонента, через скорость движения фронта кристаллизации, связан со скачком теплового потока на фазовой границе (термо - диффузионная задача Стефана) [11, т.1, с.391].

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

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

Анализ разностных методов предполагает, что по крайней мере на качественном уровне поведение решений модельной системы дифференциальных уравнений известно. Как правило, это не так. Учёт процессов на границе раздела фаз, даже в линейном приближении, не сводится к формулировке классических условий Дирихле или Неймана. Здесь возникают несамосопряжённые краевые задачи с нетрадиционными граничными условиями, свойства решений которых априори неизвестны. Поэтому исследование разностных алгоритмов в работе сочетается с анализом свойств решений модельной дифференциальной задачи. Особое внимание уделено изучению её спектра, что позволило выявить детальную картину поведения частных решений при различных значениях безразмерных параметров задачи [108], [109]. Полученные данные используются при исследовании устойчивости различных способов численной реализации граничных условий.

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

Результаты, полученные для одномерных модельных задач, позволили в классе методов с явным выделением фронта кристаллизации проанализировать поведение и область применимости основных используемых на практике алгоритмов решения задач тепло - мас-сопереноса в системах с фазовыми превращениями, указать способы отбора и построения надёжных разностных методов.

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

В термо - диффузионной задаче Стефана распределение температуры и концентрации на границе жидкое - твёрдое должны удовлетворять трём условиям [11, т.1, с.391]: два из них вытекают из законов сохранения энергии и массы, они связывают тепловые и концентрационные потоки со скоростью движения фронта; третье

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

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

В работе предложен алгоритм совместного решения уравнений переноса тепла и массы. В нём на фронте кристаллизации задаются два условия. Одно из них, по - прежнему, фазовая диаграмма, второе

- нелинейное соотношение, связывающее температуру и концентрацию. Оно является следствием двух законов сохранения на границе жидкое - твёрдое и получено путём исключения из соответствующих уравнений скорости движения границы фазового перехода. Фактически это соотношение означает, что, как и в дифференциальном случае, скорость движения границы раздела фаз, вычисленная по потоку растворённого компонента, совпадает со скоростью, найденной из условия баланса температуры. Уравнения для определения концентрации и температуры решаются совместно с помощью метода Ньютона. В рамках такого подхода удаётся сконструировать вычислительную процедуру, обладающую большим запасом устойчивости, чем раздельные алгоритмы, и гарантирующую получение надёжных результатов во всём диапазоне параметров, представляющих практический интерес.

О математическом моделировании процессов кристаллизации.

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

Тепло - массоперенос при получении полупроводниковых материалов методом Бриджмена подробно изучался в двумерном приближении Р.А.Брауном с соавторами [110] - [116]. В первых из этих работ задача рассматривалась в квазистационарной постановке. Она состоит в том, что скорость движения границы раздела фаз совпадает со скоростью протяжки ампулы, а поля температуры и скорости движения жидкости - стационарны. Моделирование проводилось для случая кристаллизации чистого вещества или расплава, легированного небольшим количеством примеси. Температура фазового перехода была постоянна и граница раздела фаз проходила по соответствующей изотерме.

В более поздних работах [114] - [115] использовалась нестационарная модель, в рамках которой изучалась кристаллизация различных полупроводниковых материалов. Концентрационная конвекция в этих работах не учитывалась и форма границы фазового перехода не связана с распределением примеси в жидкой и твёрдой фазах. Анализ условий выращивания монокристаллов (С ¿Т е) Х(Н дТ е)\-х методом Бриджмена на основе решения нестационарной задачи о фазовом переходе в бинарной системе содержится в [116], где изучено влияние тепловой и концентрационной конвекции на состав растущего кристалла, скорость роста и форму границы раздела фаз. В расчётах использовалась реальная фазовая диаграммы системы С ¿Те - НдТе.

Во всех упомянутых выше случаях вычислительный эксперимент проводился на суперкомпьюперах типа Cray с помощью различных модификаций специально сконструированного метода конечных элементов [117], [118].

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

Число работ, в которых численное исследование процессов кристаллизации проводится в нестационарном приближении, сравнительно невелико. Как правило они относятся к моделированию фазового перехода в чистых веществах или разбавленных растворах, когда поле концентрации не влияет на форму границы раздела фаз. (См., например, [124] - [126].)

Из публикаций, в которых рассматривается фазовый переход в системах с температурой кристаллизации, зависящей от состава жидкой фазы, отметим работы [127] и [128]. В первой из них продолжается исследование особенностей тепло - массопереноса при выращивании методом Бриджмена кристаллов CdZnTe большого радиуса. Авторы последовательно изучали задачу в квазистационарном приближении [129], затем в рамках нестационарного подхода без учёта [130] и с учётом влияния концентрации жидкой фазы на форму границы раздела фаз [127]. В этих работах использовался метод конечных элементов, близкий по свойствам алгоритмам [117], [118].

В [128] моделирование процесса зонной плавки проводится с помощью конечно - разностного алгоритма, представляющего собой обобщение на случай бинарного расплава метода [131] - [134]. При его построении в исходной задаче с подвижными границами делается зависящая от времени замена переменных, в результате которой расчётная область преобразуется в прямоугольник с фиксированными границами В новой системе координат разностная схема по пространству выписывается с помощью интегро - интерполяционного метода, производная по времени - сохраняется. Соответствующая система обыкновенных дифференциальных уравнений относительно всех искомых функций (температуры, вихря, функции тока и формы фронтов кристаллизации) решается методом [135], в котором порядок схемы по времени и шаг интегрирования выбирается исходя из поведения решения; итерации по нелинейности проводятся по методу Ньютона; при решении линейной системы итерируются последовательно отдельные группы уравнений. Сходимость итераций довольно медленная и с приемлемой точностью решение всей линейной системы уравнений находится примерно за 80 итераций

131].

Нестационарный метод численного исследования задач кристаллизации многокомпонентных соединений развивается в [136] - [138]. Разработаны математические модели, разностные алгоритмы и комплексы программ для моделирования процессов получения монокристаллов методами Бриджмена, движущегося нагревателя и бестигельной зонной плавки. Процесы рассматриваются в постоянном и вращающемся магнитном поле с учётом условий теплообмена излучением между ампулой и камерой - нагревателем в случае непрозрачной и полупрозрачной ампулы. Метод расчёта термо - диффузионной задачи Стефана строится в классе методов с явным выделением границы раздела фаз. Изменяющаяся во времени расчётная область отображается на фиксированный прямоугольник; схема выписывается в плоскости новых переменных на неподвижной сетке. Уравнения движения в форме "вихрь - функция тока" решаются совместно с помощью упомянутого выше (стр. 22) матричного алгоритма [103], поля температуры и концентрации определяются последовательно. В [136] при моделировании метода движущегося нагревателя уравнение для температуры решается с условием Стефана на фазовой границе; затем решается уравнения для концентрации. Её значение на границе раздела фаз вычисляется из фазовой диаграммы. Скорость движения фронта кристаллизации определяется из условия баланса массы на границе жидкое - твёрдое.

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

О численном исследовании метода жидкофазовой эпитаксии. В последние годы математическое моделирование процессов конвективного массопереноса в расплаве при получении эпитаксиальных слоёв полупроводниковых материалов всё чаще привлекает к себе внимание исследователей. Это связано не только с тем, что длительное время широко используемые одномерные диффузионные модели [139] - [149] в большинстве случаев явно не достаточны для анализа условий роста в эпитаксиальных системах, но и с разработкой новых способов эпитаксиального выращивания, при которых конвективное движение в расплаве, традиционно рассматриваемое как фактор, отрицательно влияющий на качество получаемых слоёв, используется для формирования условий роста, обеспечивающих получение материалов с заданными свойствами [150] - [152]. Кроме того, интерес к численному исследованию конвективного массопереноса в растворе - расплаве стимулируется работами, направленными на поиск средств управления процессами массопереноса в эпитаксиальных системах ( выращивание в "градиенте температуры," под действием электрического и магнитного поля, в условиях пониженной гравитации и т.д.).

Впервые, по - видимому, методами вычислительного эксперимента конвективный массоперенос в расплаве при получении двойных и тройных полупроводниковых соединений изучался в работах автора [98] - [102], [154]. Моделирование осуществлялось на основе решения двумерных нестационарных уравнений концентрационной конвекции, при этом использовались построенные автором данной работы численные методы.

Для систем с вертикальным расположением подложек при различных значениях параметров технологического процесса (скорости охлаждения расплава, величины зазора между подложками и др.) изучалось влияние конвективного массопереноса на скорость роста и однородность по толщине эпитаксиальных слоёв фосфида галлия. Рассматривался также процесс травления подложек в результате контакта с перегретым относительно температуры ликвидуса раствором [98], [100], [102].

На основе анализа результатов расчётов и экспериментальных исследований, проведённых Э.А.Твировой и Л.А.Дмитриевой в институте ГИРЕДМЕТ, были предложены технологические режимы, обеспечивающие получение однородных по толщине эпитаксиальных слоёв при высоких скоростях охлаждения и большой величине зазора между подложками. Найденный таким образом способ выращивания был зарегистрирован как изобретение [153].

В горизонтальных системах изучался процесс конвективного массопереноса в расплаве при выращивании слоёв арсенида галлия на две подложки, расположенные над расплавом и под ним [99], [102]. В зависимости от величины зазора и длины подложек в расчётах были получены различные режимы течения и изучена их связь с формой поверхности ЭС.

Численному исследованию эпитаксиального выращивания тройных полупроводниковых соединений посвящена работа [154], где рассматривается процесс получения гетероструктур С¿¿уНС\уТе. В работе анализуется влияние различных параметров технологического процесса на динамику роста, геометрические параметры и состав эпитаксиальных слоёв, однородность распределения состава по толщине и площади монокристалла; рассмотрирается важный с практической точки зрения вопрос о подрастворении подложки на начальных стадиях процесса, указаны технологические режимы, обеспечивающие получение материалов с заданными свойствами.

Математическое моделирование конвективного массопереноса при получении двойных соединений методом ЖФЭ проводится также в [155] - [158], где значительное внимание уделяется параметрическому анализу процесса и поиску режимов, способствующих уменьшению геометрической неоднородности эпитаксиальных структур, в частности моделируется рост слоёв в невесомости, при медленном вращении подложек и т.д.

Моделированию процесса ЖФЭ посвящён недавний цикл работ [159] - [165], где вычислительный эксперимент последовательно применяется для анализа условий роста двойных и тройных полупроводниковых соединений. Задача решается методом конечных разностей или методом конечных элементов [160]; подробно обсуждаются трудности, связанные с численной реализацией условий на границе раздела фаз [164].

Содержание диссертации по главам.