автореферат диссертации по информатике, вычислительной технике и управлению, 05.13.18, диссертация на тему:Численное исследование конвективной устойчивости при выращивании эпитаксиальных слоев

Введение.

1 Алгоритм численного решения трехмерных уравнений Навье-Стокса для несжимаемой жидкости в естественных переменных

1.1 Система уравнений.

1.2 Расщепление по физическим процессам.

1.3 Численный метод решения уравнений Навье-Стокса для несжимаемой жидкости.

1.4 Анализ различных способов решения уравнения для давления.

1.5 Тестовые расчеты

1.6 Основные результаты.

2 Конвективная устойчивость

2.1 Основные предположения.

2.2 Математическая модель.

2.3 Метод решения задачи концентрационной конвекции

2.4 Модельная задача.

2.5 Некоторые дополнительные теоретические и экспериментальные данные

2.6 Численное исследование возникновения подкритической конвекции.

2.7 Зависимость структуры течения от величины зазора

2.8 Влияние конвекции на форму поверхности эпитакси-ального слоя

2.9 Сравнение с результатами двумерных расчетов.

2.10 Анализ влияния параметров сетки на результаты расчетов

2.11 Основные результаты

3 Математическое моделирование технологических режимов выращивания ЭС с предварительным подрастворе-нием подложки

3.1 Описание процесса.

3.2 Результаты расчетов.

3.3 Учет конечной теплопроводности стенок ростовой камеры и жидкости.

3.4 Основные результаты

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

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

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

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

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

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

Численное моделирование процессов конвективной диффузии в растворе-расплаве при жидкофазовой эпитаксии ранее основывалось преимущественно на двухмерных моделях. Целью данной работы является изучение в трехмерной постановке гидродинамических процессов протекающих в расплаве при ЖФЭ, определение порога устойчивости системы, динамики развития неустойчивостей, анализ влияния возникающего конвективного движения на характеристики ЭС. Также исследуется зависимость наблюдаемых в расчетах течений и характеристик эпитаксиальных слоев от неоднородности распределения температур в ростовой камере. Проводится сравнение результатов, полученных в рамках двумерного и трехмерного приближений, определяется диапазон параметров, в котором существенно использование трехмерных моделей.

Процесс развития конвективной неустойчивости рассматривается в приближении Обербека-Буссинеска [2, 3]. Соответствующая система уравнений получается из уравнений Навье-Стокса в предположении, что плотность жидкости не зависит от давления и является линейной функцией температуры и концентрации. Изменение плотности учитывается только в уравнении импульса в слагаемом, связанном с подъемной силой; градиенты скорости достаточно малы и тепловыделением за счет вязкой диссипации можно пренебречь. В результате процесс описывается системой уравнений Навье-Стокса для несжимаемой жидкости с подъемной силой в уравнении для импульса, и уравнениями конвективной теплопроводности и диффузии. Важную роль при моделировании процессов конвективного тепломассопереноса играет выбор алгоритма решения уравнений Навье-Стокса для несжимаемой жидкости.

О методах численного решения уравнений Навье-Стокса для несжимаемой жидкости.

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

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

Для численного моделирования несжимаемых течений применяются два подхода. Первый из них состоит в том, что из исходной системы уравнений исключается давление и задача записывается в иных переменных: "функция тока-вихрь" [4,5], "вихрь-скорость" [4] или "функция тока-скорость" [6]. Второй подход основан на использовании естественных переменных скорость-давление [4,7-9].

Подход, основанный на использовании переменных "функция тока-вихрь", активно применяется для решения двумерных задач. Его достоинством является автоматическое выполнение условия несжимаемости, меньшее число неизвестных и отсутствие сложностей, связанных с определением давления. При этом, однако, возникают нелокальные граничные условия для вихря, строгая реализация которых требует совместного решения уравнений для вихря и функции тока или скорости [10-14]. Обобщение этого подхода на трехмерный случай сопряжено с немалыми трудностями и проигрывает алгоритмам в естественных переменных как по удобству использования, так и по числу уравнений [15]. Поэтому в трехмерном случае основные усилия направлены на построение численных методов в переменных скорость-давление.

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

Основным требованием, предъявляемым к дискретной модели, является выполнение законов сохранения, справедливых в исходной физико-математической постановке задачи, так называемый принцип консервативности [16-18]. Уравнения Навье-Стокса для несжимаемой жидкости представляют собой записанные в дифференциальной форме законы сохранения массы (уравнение неразрывности) и импульса. Уравнение для кинетической энергии отсутствует, а баланс кинетической энергии является следствием выполнения законов сохранения массы и импульса. В разностном случае выполнение законов сохранения массы и импульса достигается аппроксимацией дивергентной формы исходных дифференциальных уравнений интегро-интерполяционным методом [17]. Однако численный метод должен сохранять не только массу и импульс, но и кинетическую энергию. Это особенно важно при моделировании существенно нестационарных задач на больших промежутках времени. В разностном случае из выполнения балансов массы и импульса, вообще говоря, не следует выполнения баланса кинетической энергии. Можно построить схемы, удовлетворяющие законам сохранения массы и импульса, но не сохраняющие кинетическую энергию [19,20].

Применительно к уравнениям Навье-Стокса для несжимаемой жидкости необходимость выполнения принципа консервативности впервые продемонстрировал А. Агака\уа на примере алгоритмов в переменных "функция тока-вихрь" [21]. При интегрировании уравнений на больших временных отрезках нарушение в дискретной модели баланса кинетической энергии приводило к ее неограниченному росту. Аналогичное явление наблюдал Б.Л. Рождественский, используя гибрид метода конечных разностей и спектрального метода для моделирования турбулентности [22].

Если в дискретной модели закон сохранения энергии не выполнен, в системе наблюдается либо избыточная диссипация, либо нарастание кинетической энергии. В случае избыточной диссипации схема устойчива, но результаты, полученные с ее помощью при высоких значениях числа Рейнольдса, нефизичны: наблюдается сильное затухание энергии в областях высоких градиентов. Такими, например, являются схемы с аппроксимацией конвективных членов против потока первого порядка [19]. Нарастание кинетической энергии приводит к неустойчивости численной процедуры. Особенно опасны схемы, устойчивые при небольших значениях числа Рейнольдса и неустойчивые при больших. В расчетах с умеренным значением числа Рейнольдса они могут оставаться устойчивыми, но давать существенную погрешность за счет неконсервативности [19,23].

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

Сравнение алгоритмов на разнесенных и совмещенных сетках, выполненное в [26] для случая сжимаемых течений, показало, что на разнесенных сетках лучше выполняются законы сохранения. Другим преимуществом схем на разнесенных сетках является связь между скоростью и давлением, не дающая нефизических, сеточных осцилляций, характерных для схем на совмещенных сетках [20,27].

Схемы порядка точности выше второго и сохраняющие кинетическую энергию, можно найти в работах [19] и [20] для случая разнесенных и совмещенных сеток соответсвенно. Также были разработаны компактные схемы высокого порядка точности [28]. S.K. Lele показал, что компактные схемы высокого порядка точности имеют меньший размер шаблона, и при определенном выборе параметров, лучше разрешают высокочастотные гармоники, чем стандартные конечно-разностные схемы того же порядка точности. Первоначально компактные схемы строились на основе формальной аппроксимации дифференциальных уравнений. В более поздних работах был применен интегро-интерполяционный метод [29]. Схемы высокого порядка точности, благодаря малой схемной вязкости, подвержены осцилляциям, приводящим к вычислительной неустойчивости [26,28,30] и, чем менее подробная сетка используется, тем выше вероятность их возникновения. Для борьбы с осцилляциями приходится вводить искусственную вязкость, фильтрующую высокочастотные гармоники, как, например, предложено в [28], или использовать регуляри-заторы высокого порядка точности [24,31-33].

Большинство методов решения уравнений Навье-Стокса для несжимаемой жидкости основано на дискретизации по времени типа предиктор-корректор [17,34,35]. Впервые алгоритмы такого типа были использованы в работах [36-38]. Их суть состоит в том, что на каждом временном слое сначала вычисляется предиктор скорости, вообще говоря, не удовлетворяющий условию несжимаемости. Затем из уравнения Пуассона определяется давление. При помощи найденного давления предиктор скорости корректируется таким образом, чтобы удовлетворить уравнению неразрывности. Это можно интерпретировать как проекцию предиктора скорости на множество соленоидальных функций, поэтому такие методы называют еще проекционными [37, 38]. Полученное после этапа коррекции поле скоростей считается найденным на текущем временном слое и осуществляется переход к следующему шагу по времени. В [36] использовалась разнесенная сетка, в [37,38] неразнесенная. Основные черты этих методов сохраняются до сих пор. Изложенные алгоритмы могут рассматриваться также как частный случай метода дробных шагов Яненко [34].

В проекционных методах на промежуточных этапах вычислений требуются краевые условия для предиктора скорости, в свою очередь определяющие краевые условия для давления. В 1972 С. Hirt и J. Cook [39] предложили использовать метод [36] с уравнением Пуассона для определения поправки давления на следующем временном слое, вместо вычисления самого давления. Это позволяет достаточно просто ставить краевые условия для предиктора скорости со вторым порядком точности по времени [23]. Проблема выбора краевых условий обсуждается в [4,40], где показано, что краевые условия для давления, отсутствующие в дифференциальной постановке задачи, порождают ошибку, имеющую характер пограничного слоя, толщина которого пропорциональна корню из произведения вязкости на величину шага по времени [41].

В работах [41-45] показано, что порядок точности вычисления скорости на следующем временном слое совпадает с порядком аппроксимации метода расщепления. Точность определения давления при этом имеет лишь первый порядок по времени [46]. Использование явных формул, применявшихся в ранних методах в расчете предиктора скорости, приводит к серьезным ограничениям на величину шага по времени: 0.25У2гЯе < 1 и г/(Яе/г2) < 0.25 [4,47,48], избежать которых удается используя неявные и полунеявные схемы.

Существуют также алгоритмы, основанные на интегрировании уравнений Навье-Стокса по времени методом Рунге-Кутты [49] и Рунге-Кутты с совместным решением получающихся нелинейных уравнений итерациями по Ньютону [50]. Достоинством такого подхода является отсутствие необходимости ставить краевые условия для давления, что полностью согласуется с физической постановкой задачи. Однако решение нелинейных уравнений более трудоемко и сходимость таких методов часто вызывает затруднения. Попытка построить конечно-разностный метод, в которым система уравнений решается относительно вектора неизвестных, содержащего все компоненты скорости и давление, была предпринята в [48]. При решении двумерных задач, алгоритм оказался безусловно устойчивым, однако его обобщение на трехмерный случай связано со значительными трудностями, возникающими при решении соответствующей системы сеточных уравнений.

Для ускорения расчетов, ведутся разработки параллельных алгоритмов численного исследования процессов описываемых нестационарными уравнениями Навье-Стокса [51].

В данной работе для решения уравнений Навье-Стокса используется вычислительная процедура типа предиктор-корректор [7] второго порядка точности по пространству и первого по времени. Пространственная аппроксимация не вносит вклада в баланс кинетической энергии [24,52]. Метод использует центральные разности на разнесенной сетке. Давление относится к центрам ячеек, скорости - к центрам граней [53].

Неустойчивость Рэлея-Бенара

Начало систематическому исследованию конвективной неустойчивости положили экспериментальные работы Бенара. Он наблюдал возникновение движения в тонком, горизонтальном, подогреваемом снизу слое жидкости со свободной поверхностью [54]. Первые попытки теоретического исследования такой задачи были предприняты лордом Рэлеем, который провел линейный анализ задачи с двумя свободными границами [55]. Поэтому движение, возникающее в тонком горизонтальном слое жидкости под действием градиента температуры или концентрации, получило название конвекции Рэлея-Бенара. В настоящее время, этому вопросу посвящено огромное число работ, однако теория конвективной неустойчивости еще далека от своего завершения, ограничиваясь анализом частных задач и сопоставлением полученных результатов с экспериментальными данными [3,56-58].

Неустойчивость Рэлея-Бенара является одной из центральных проблем механики жидкости. Интерес к этому явлению объясняется тем, что в рамках простой физической постановки задачи удается получить широкий класс течений: от устойчивых и регулярных, до нерегулярных, быстро изменяющихся во времени. При этом в широком диапазоне параметров в толщине слоя сохраняется простая структура. Основные вопросы, возникающие при исследовании конвекции Рэлея-Бенара, это: изучение порога устойчивости рассматриваемой системы и планформы возникающего движения (проекции конвективной структуры на горизонтальную плоскость); анализ влияния структуры течения на характеристики тепломассопереноса; определение скорости развития конвективного движения из состояния покоя [3,58].

Экспериментальное исследование указанных явлений сопряжено с определенными трудностями. В первую очередь они связаны с большими временами эволюции конвективных структур (для воды от 10 минут до сотен часов, в зависимости от горизонтальных размеров области), высокой стоимостью экспериментальных установок, особенно для больших аспектных отношений, необходимых для полноценного анализа конвективных структур [59]. Шагом вперед в этой области было применение жидких газов, позволивших снизить характерные времена и расширить диапазон изменения значения числа Прандля [59]. Также представляет сложность само наблюдение за экспериментом и анализ структуры движения по полученным измерениям. Поэтому проведение численного моделирования и последующее сравнение полученных результатов с контрольными экспериментами является эффективным способом изучения конвективной неустойчивости.

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

Найти и исследовать стационарные течения в надкритической области позволяет слабо-нелинейный анализ, основанный на представлении решения в виде ряда по степеням малого амплитудного параметра [3,63,64]. Этот параметр пропорционален корню квадратному из над-критичности: разности между значением числа Рэлея и его критической величиной. Сам по себе расчет стационарных надкритических движений не дает возможности ответить на вопрос, какое именно движение реализуется в действительности. В экспериментах могут наблюдаться лишь устойчивые, по меньшей мере к бесконечно-малым возмущениям, течения. Поэтому найденные решения исследуются на устойчивость при помощи линейного анализа. Наиболее полные результаты по исследованию устойчивости надкритических течений были получены А. Шлютером, Д. Лорцем и Ф. Буссе. Они показали, что в областях небольшой надкри-тичности устойчивыми относительно трехмерных возмущений являются только двумерные валы [64,65]. В дальнейшем были найдены диаграммы устойчивости различных типов течений, получившие название баллонов или сачка Буссе [62,65,66]. Однако механизм отбора планформ и волновых чисел до сих пор остается неопределенным [58].

Анализ влияния боковых стенок на структуру течения позволяет сузить спектр допустимых течений, сделав его дискретным. Но и в этом случае вопросы устойчивости и реализуемости течений остаются открытыми. При достаточно большом аспектном отношении влияние боковых стенок на критическое значение числа Рэлея мало, основное влияние стенок на структуру течения наблюдается в непосредственной близости от них [67-69]. Течение, как правило, имеет двумерную структуру, валы выстраиваются параллельно короткой стороне области или подходят к границам практически под прямым углом [68-70].

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

Более важным для рассматриваемого в данной работе класса задач является случай, когда верхняя и нижняя границы области равномерно охлаждаются с постоянной скоростью. В результате внутри области формируется нелинейный профиль температуры. Эта задача эквивалентна стационарной задаче об устойчивости неподвижной жидкости в бесконечном горизонтальном слое с параболическим профилем температуры [71]. Нелинейный анализ, основанный на разложении в ряд по степеням малых параметров: амплитудного и кривизны профиля температуры, показал, что критическое значение числа Рэлея в данном случае ниже, чем при постоянном градиенте, и зависит от кривизны профиля [71]. Как и при переменной вязкости, здесь возможно возникновение конечно-амплитудной подкритической конвекции в форме шестиугольных ячеек, направление движения жидкости в которых определяется знаком кривизны профиля температуры. Эти задачи объединяет отсутствие симметрии верх-низ; для них характерны шестиугольные ячейки, в то время как для задач симметричных относительно средней плоскости, например при линейном профиле, типичны валы.

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

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

Численное моделирование конвекции в трехмерной постановке началось достаточно давно. В [72] исследуется процесс формирования конвективных структур на основе двумерной модели трехмерной конвекции; описываются процессы возникновения и взаимодействия конвективных структур, а также их дефектов. В [73] для областей с небольшим аспектным отношением (<3) в приближении бесконечного числа Прандтля изучается устойчивость двумерных валов относительно трехмерных возмущений. Численное исследование устойчивости жидкости с вязкостью, зависящей от температуры, для бесконечного значения числа Прандтля проведено в [74], где, в частности, было получено движение, охватывающее только часть слоя, с планформой в виде квадратных ячеек. Важность учета трехмерности показана в [15] на примере задаче о тепловой гравитационной конвекции в области с боковым подогревом, и в [75,76]. С увеличением значения числа Рэлея течение теряет двухмерную структуру и становится трехмерным. В [77] проводится численное моделирование подкритической конвекции в жидкости с вязкостью, зависящей от температуры. Для этого в шестиугольной планформе выбирается одна ячейка, ее внутренняя часть приближается цилиндром. На границе цилиндра ставятся периодические краевые условия и, используя предположение осесимметричности, задача сводится к двумерной. В работе получена зависимость критического числа Рэлея от перепада вязкости в области. В последнее время ведется интенсивное исследование конвекции Рэлея-Бенара в вязком сжимаемом газе [78-80]. Также, созданы пакеты программ, позволяющие осуществлять моделирование процессов тепломассообмена [81].

В данной работе изучается процесс потери конвективной устойчивости раствора-расплава при выращивании трехкомпонентных монокристаллических слоев методом ЖФЭ. Процесс возникновения и развития конвекции в системе во многом аналогичен случаю, рассмотренному в [71], но существенно нестационарен, т.к. граничные условия зависят от времени. На основе анализа результатов двумерных расчетов найдено значение критического числа Рэлея. В трехмерных расчетах зарегистрировано существование подкритической конвекции в форме шестиугольных ячеек; определено минимальное значение числа Рэлея, при котором наблюдается конвективное движение. Полученные результаты хорошо согласуются с теоретическими данными [71].

Учет теплопроводности стенок.

Влияние конечной теплопроводности стенок на величину критического числа Рэлея было отмечено еще в 1926 году Н. Jeffreys [82]. В 1964 Е.М. Sparrow с соавторами выполнили анализ процесса развития неустойчивости в бесконечном горизонтальном слое, теплопередача на верхней и нижней границе которого определялась линейным законом Фурье, т.е. краевыми условиями третьего рода. Для нескольких случаев было определено критическое значение числа Рэлея в виде функции числа Био (характеризующего отношение теплопроводности внутри области и во внешней стреде) [83]. D.T. Hurle и Г.З. Гершуни с соавторами изучали устойчивость равновесия бесконечного горизонтального слоя жидкости, ограниченного сверху и снизу полубесконечными стенками, теплопроводность которых различна и отличается от теплопроводности жидкости [3,84]. Было показано, что чем ниже теплопроводность стенок, тем ниже порог устойчивости: он опускается с 1708 для идеально проводящих стенок до 720 для теплоизолированных. Волновое число при этом уменьшается с 3.11 до нуля.

Исследованию влияния боковых стенок на порог устойчивости посвящены работы [67,85-88]. Основные результаты этих исследований состоят в следующем: теплопроводность боковых стенок влияет на значение критического числа Рэлея только при аспектном отношении меньше двух, при этом значение критического числа Рэлея для области с непроводящими стенками значительно ниже, чем в случае идеально теплопро-водящих боковых границ. Очевидно, что первый результат обусловлен вязким взаимодействием жидкости с боковыми стенками, а второй - тем что температурные флуктуации вблизи идеально теплопроводной стенки затухают быстрее, чем возле теплоизолированной. Эксперименты, проведенные K.Stork и U.Müller'oM для области с различной теплопроводностью верхней и нижней стенок и с плохо проводящими боковыми границами подтвердили эти результаты [68].

Работы по изучению влияния конечной теплопроводности боковых стенок на значение критического числа Рэлея и реализующиеся формы течений продолжаются и в настоящее время, однако вопрос в целом все еще остается открытым. Из более поздних работ отметим [89], в которой рассматривается бесконечный горизонтальный слой жидкости, ограниченный стенками конечной толщины и теплопроводности. Расчет линеаризованной задачи показал, что при отношении толщины горизонтальных границ к толщине жидкой фазы больше единицы, результаты близки к полученным Г.З. Гершуни для полубесконечных стенок [3]. Если же это отношение порядка единицы, то для стенок с теплопроводностью на порядок выше чем в жидкости, разница между двумя результатами достигает 20%. Для волнового числа это влияние более существенно. При отношения толщин порядка единицы, разница между результатами [89] и [3] в зависимости от теплопроводностей стенок и жидкости может достигать 20%, если стенки в 10 раз тоньше слоя жидкости - то разница может возрастать до 100%. В [89] также отмечалось, что даже небольшие перепады температуры у границы, являются причиной возникновения пристеночного вихря.

Таким образом, учет конечной теплопроводности и толщины стенок может оказать существенное влияние на наблюдаемые картины течений в задаче Рэлея-Бенара. В данной работе проводится численное исследование процесса роста кристалла методом ЖФЭ с учетом конечной теплопроводности стенок ростовой камеры. Результаты расчетов сравниваются с полученными для случая идеально-проводящих стенок.

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

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

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

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

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

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

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

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

Далее, проводится численное моделирование процесса роста ЭС в трехмерной постановке для значений числа Рэлея в диапазоне 1.1-103 — 4.5-105. Изменение значения числа Рэлея осуществлялось за счет выбора толщины слоя жидкости. При моделировании обнаружена подкрити-ческая конвекция, которая является существенно трехмерной, и имеет форму шестиугольных ячеек. Форма подкритического движения и минимальное значение числа Рэлея, при котором оно наблюдается, полностью согласуются с результатами конечно-амплитудной теории устойчивости [3,58,71]. В надкритической области, для значений числа Рэлея в диапазоне 1.6Яа£Т < 11а < 5Яасг, Яа^ = 1545 зарегистрирована устойчивая ячейковая конвекция, при 5Ла^. < Ла < 2(Шасг спицевидная, и для Ла ^ 20Ласг хаотическая.

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

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

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

• в областях большой надкритичности, скорость перестройки трехмерных течений выше, а неоднородность поверхности выросшего слоя - меньше, чем в двумерном.

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

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

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

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

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

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

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

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

Основные результаты работы:

1. Разработан разностный алгоритм и комплекс программ для численного исследования трехмерных нестационарных процессов конвективного тепломассопереноса применительно к проблеме выращивания трехкомпонентных эпитаксиальных слоев из жидкой фазы.

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

3. В диапазоне изменения параметров, представляющем практический интерес, найдены технологические режимы, при которых мас-соперенос в жидкой фазе осуществляется преимущественно механизмом диффузии. В системах с предварительным растворением подложки проанализировано влияние динамики развития возмущений полей температуры и концентрации на характеристики ЭС. Показано, что скорость роста амплитуды возмущений в форме двумерных валов выше, чем у трехмерных ячеек, и они приводят к большей неоднородности поверхности слоя. Перепад температуры порядка 0.1-0.3 градуса при абсолютной величине Т ~ 700.fi, порождаемый конечной теплопроводностью веществ, приводит к увеличению неоднородности толщины выросшего слоя по сравнению с идеально теплопроводящей средой в 1.5-2 раза.

Автор выражает искреннюю благодарность: своему научному руководителю д.ф.-м.н. Мажоровой О.С. и член, корреспонденту РАН, д.ф.-м.н. Попову Ю.П. за постоянное внимание и помощь в работе, участие в многочисленных обсуждениях; сотрудникам 27-й лаборатории института ГИРЕДМЕТ за постановку задачи; всем сотрудникам 11-го отдела ИПМ им М.В. Келдыша РАН за сотрудничество и ценные советы во время подготовки диссертационной работы.

3.4 Основные результаты

1. Изучено влияние конвективного движения на процесс растворения подложки в методе жидкофазовой эпитаксии. Показано, что начальные возмущения в форме двумерных валиковых структур обладают большей скоростью роста, чем квадратные ячейки. При малой надкритичности (Яа ~ 2.7Яасг), это приводит к большей неоднородности поверхности травления. Поэтому они являются наиболее неблагоприятными с технологической точки зрения. При средней надкритичности (Яа ~ 4.3Дасг) неоднородность поверхности фронта травления практически не зависит от формы вносимого начального возмущения. Это позволяет при умеренном значении числа Яа < 5Яасг использовать двумерные модели для анализа задачи растворения подложки.

2. Наблюдаемые в расчетах существенно трехмерные структуры течения образованы ячейками g-типa, что полностью согласуется с данными теоретического анализа [71].

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

1. Spiegel Е.А., Veronis G. On the Boussinesq approximation for a compressible fluid. Astrophys. J., 131(2), 442-447, 1960.

2. Гершуни Г.З., Жуховицкий E.M. Конвективная устойчивость несжимаемой жидкости. М.: Наука, 1972. 392 стр.

3. Флетчер К. Вычислительные методы в динамике жидкости. М.: Мир, 1991. т. 1. 502 е., т.2. 552 с.

4. Пасконов В.М., Полежаев В.П., Чудов JI.A. Численное моделирование процессов тепло- и массообмена. М.: Наука, 1984. 288 стр.

5. Gupta М.М., Kalita J.С. A new paradigm for solving Navier-Stokes equations: streamline-velocity formulation. J. Comput. Phys., 207,5268, 2005.

6. Белоцерковский O.M. Численное моделирование в механике сплошных сред. М.: ФизМатЛит, 1994.

7. Patankar S.V. Numerical heat transfer and fluid flow. Hemisphere Publishing Corporation, 1996.

8. Ferziger J.H., Peric M. Computational Methods for Fluid Dynamics. 2002.

9. Quartapelle L., Valz Gris F. Projection condition on the vorticity in viscous incompressible flows. Int. J. for numerical methods in fluids, 1, 129-144, 1981.

10. Мажорова О.С., Попов Ю.П. Об одном алгоритме численного решения двумерных уравнений Навъе Стожа, препринт Института прикладной математики им. М.В.Келдыша РАН, - Москва., 1979, № 37. 24 с.

11. Мажорова О.С., Попов Ю.П. О методах численного решения уравнений Навье Стокса. Ж. Вычисл. Матем и мат. физики, 20(4), 1005-1020, 1980.

12. Мажорова О.С., Попов Ю.П. Матричный итерационный метод численного решения двумерных уравнений Навье Стокса. ДАН СССР, 259(3), 535-540, 1981.

13. Бессонов O.A., Брайловская В.А., Полежаев В.И. Пространственные эффекты конвекции в расплавах: концентрационные неоднородности, возникновение несимметрии и колебания. Изв. РАН МЖГ., 3, 74-82, 1997.

14. Тихонов А.Н., Самарский A.A. О сходимости разностных схем в классе разрывных коэффициентов. ДАН СССР, 124(3), 1959.

15. Самарский A.A. Теория разностных схем. М.: Наука, 1989.

16. Самарский A.A., Попов Ю.П. Разностные методы решения задач газовой динамики. Изд.4- УРСС, 2004.

17. Morinishi Y., Lund T.S., Vasilyev O.V., Moin P. Fully conservative higher order finite difference schemes for incompressible flow. J. Comput. Phys., 143, 90-124, 1998.

18. Vasilyev O.V. High-order finite difference schemes on non-uniform meshes with good conservation properties. J. Comput. Phys., 157, 746-761, 2000.

19. Arakawa A. Computational design for long-term numerical integration of the equation of fluid motion: Two dimentional incompressible flow. J. Comput. Phys., 1, 119-143, 1966.

20. Рождественский Б.JI., Моисеенко Б.Д., Сидорова В.К. Условия численного моделирования предельных режимов течений вязкой жидкости, препринт Института прикладной математики им. М.В.Келдыша РАН, Москва., 1979, № 14.

21. Kim J., Moin P. Application of a Fractional-Step Method to Incompressible Navier-Stokes Equations. J. Comput. Phys., 59, 308323, 1985.

22. Гончаров A.JI., Фрязинов И.В. Разностные схемы на девтиточеч-ном шаблоне „крест"для решения уравнений Навье-Стокса. Ж. Вычисл. Матем и мат. физики, 28, 867-878, 1988.

23. Вабищевич П.Н., Павлов А.Н.,Чурбанов А.Г. Методы расчета нестационарных несжимаемых течений в естественных переменных на неразнесенных сетках. Мат. Модел., 8(7), 81-108, 1996.

24. Nagarajan S., Lele S.K., Ferziger J.H. A robust high-order compact method for large-eddy simulation. J. Comput. Phys., 191, 392-419, 2003.

25. Gresho M.P., Sani R.L. On pressure boundary conditions for the incompressible Navier-Stokes equations. Int. J. Numer. Methods. Fluids, 7, 1111-1145, 1987.

26. Lele S.K. Compact finite difference schemes with spectral-like resolution. J. Comput. Phys., 103, 16-42, 1992.

27. Kobayashi M.H. On a class of Pade finite volume methods. J. Comput. Phys., 156, 137-180, 1999.

28. Piller M., Stalio E. Finite-volume compact schemes on staggered grids. J. Comput. Phys., 197, 299-340, 2004.

29. Гончаров A.Jl., Фрязинов И.В. О построении монотонных разностных схем для уравнений Навъе Стокса на девятиточечных шаблонах, препринт Института прикладной математики им. М.В.Келдыша РАН, - Москва., 1986, № 93. 14 с.

30. Гончаров A.JL, Фрязинов И.В. Разностные схемы на девятиточечных шаблонах „крест" для уравнений Навъе Стокса в переменных скорость - давление, препринт Института прикладной математики им. М.В.Келдыша РАН, - Москва., 1986, N2 53. 17 с.

31. Мажорова О.С., Марченко М.П., Фрязинов И.В. Монотонизи-руещее регуляризаторы и матричный метод решения уравнений Навье-Стокса для несжимаемой жидкости. Мат. Модел., 6(12), 97-116, 1994.

32. Ковеня В.М., Яненко Н.Н. Методы расщепления в задачах газовой динамики. Н.: Наука, 1981.

33. Марчук Г.И. Методы расщепления. М.: Наука, 1988.

34. Harlow F.H., Welch J.E. Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface. Physics of Fluids, 8(12), 2182-2189, 1965.

35. Chorin A.J. Numerical solution of Navier-Stokes equations. Math.Comput., 22, 745-762, 1968.

36. Temam R. Une méthode d'approximations de la solution des équations de Navier-Stokes. Bull. Soc. Math. France 98, стр. 115-152, 1968.

37. Hirt C.W., Cook J.L. The Calculation of Three-Dimensional Flows Around Structures and Over Rough Terrain. J. Comput. Phys., 10, 324-340, 1972.

38. Orszag S.A., Israeli M., Deville M.O. Boundary conditions for incompressible flows. J. Sci. Comput., 1, 75-111, 1986.

39. Blasco J., Codina R. Error estimates for an operator-splitting method for incompressible flows. Applied Numerical Mathematics, 51, 1-17, 2004.

40. Shen J. On error estimates of projection methods for Navier-Stokes equations: First-order schemes. SIAM J Numer Anal, 29, 57-77, 1992.

41. Shen J. On error estimates of projection methods for Navier-Stokes equations: second-order schemes. Math. Comput., 65, 1039-1065, 1996.

42. Guermond J.L., Quartapelle L. On the approximation of the unsteady Navier-Stokes equations by finite elements projection methods. Numer. Math., 80, 207-238, 1998.

43. Ying Lomg-an. Viscosity-splitting scheme for the Navier-Stokes equations. Numer methods Partial Differential Equations, 7, 317-338, 1991.

44. Perron S., Boivin S., Herard J-M. A finite volume method to solve the 3D Navier-Stokes equations on instructured coolocated meshes. Int. J. Numer. Methods. Fluids, 33, 1305-1333, 2004.

45. Peyret R., Taylor T.D. Computational Methods for Fluid Flow. 1983.

46. Курганов Д.В., Мажорова О.С., Попов Ю.П. О методах решения уравнений Навье-Стокса в естественных переменных, препринт Института прикладной математики им. М.В.Келдыша РАН, Москва., 2001, № 39. 31 с.

47. Nikitin N.V. Third-order-accurate semi-implicit Runge-Kutta scheme for incompressible Navier-Stokes equation. Int. J. Num. Meth. Fluids, 51(2), 2005.

48. Pereira J.M.C., Kobayashi M.P., Pereira J.C.F. A fourth-order accurate finite volume compact method for incompressible Navier-Stokes solutions. J. Comput. Phys., 167, 217-243, 2001.

49. Горбунов А.А, Никитин С.А., Полежаев В.И. Численное моделирование конвекции Рэлея-Бенара в сжимаемом газе на параллельных вычислительных системах. Вестник ННГУ. Математическое моделирование и оптимальное управление, 1(28), 75-82, 2005.

50. Гончаров А.Л., Фрязинов И.В. Сеточный метод решения трёхмерных уравнений Навье Стокса в параллелепипеде. Дифференц. уравн., 27(7), 1137-1144, 1991.

51. Колмычков В.В., Мажорова О.С., Попов Ю.П. К расчету уравнений Навье-Стокса в естественных переменных, препринт Института прикладной математики им. М.В.Келдыша РАН, Москва., 2001, № 60, 39 стр.54