Содержание к диссертации
Введение
1. Обзор литературы по изучаемой проблеме 12
1.1. Выращивание полупроводниковых кристаллов из расплава по методу Чохральского: введение в проблему 12
1.2. Пространственный характер конвекции расплава при выращивании кристаллов по методу Чохральского 17
1.2.1. Экспериментальные данные 17
1.2.2. Результаты численного моделирования 27
1.3. Основные механизмы неустойчивости при конвекции расплава 35
1.3.1. Неустойчивость подогреваемого снизу вращающегося слоя 36
1.3.2. Проявление бароклинной неустойчивости 39
1.4. Выводы 43
2. Математическая модель и численный метод 45
2.1. Определяющие уравнения 45
2.2. Численный метод и его реализация 48
2.2.1. Общая характеристика численного метода 48
2.2.2. Преобразование координат 50
2.2.3. Геометрия ячеек 54
2.2.4. Дискретизация определяющих уравнений 57
2.2.4.1. Дискретизация по времени и расчёт поправок 57
2.2.4.2. Пространственная дискретизация и расчёт невязок 63
2.2.5 Использование блочно-структурированных сеток 68
3. Методические расчёты 70
3.1. Конвекция в квадратной полости 70
3.1.1. Предварительные замечания 70
3.1.2. Постановка задачи 73
3.1.3. Стационарная конвекция при Рг= 0,71 74
3.1.4. Нестационарная конвекция приРг = 0,71 77
3.1.5. Стационарная конвекция приРг = 0,023 82
3.1.6. Нестационарная конвекция при Pr = 0,023 83
3.2. Конвекция в кубической полости 91
3.2.1. Предварительные замечания 91
3.2.2. Постановка задачи 92
3.2.2. Результаты расчётов конвекции при Рг = 0,71 93
3.3. Конвекция в модельной установке метода Чохральского 96
3.3.1. Предварительные замечания и постановка задачи 96
3.3.2. Нестационарная конвекция при малой надкритичности 97
4. Конвекция в емкостях упрощённой геометрии 106
4.1. Проявление бароклинной неустойчивости
во вращающихся кольцевых полостях 106
4.1.1. Сопоставление с экспериментами 107
4.1.1.1. Постановка задачи 107
4.1.1.2. Регулярные волны в воде 108
4.1.1.3. Негеострофическая турбулентность при течении ртути 117
4.1.2. Конвекция расплава кремния 124
4.1.2.1. Постановка задачи 124
4.1.2.2. Негеострофическая турбулентность в расплаве кремния 125
4.2. Проявление бароклинной неустойчивости во вращающейся цилиндрической ёмкости 133
4.2.1. Предварительные замечания и постановка задачи 133
4.2.2. Развитие регулярной волновой структуры 135
4.3. Неустойчивость центрального вихря во вращающихся цилиндрических емкостях: модельная задача в отсутствие сил плавучести 145
4.3.1. Предварительные замечания и постановка задачи 145
4.3.2. Развитие неустойчивости центрального вихря: квазипериодический и стохастический режимы 147
5. Конвекция в емкостях с геометрией, типичной для тиглей метода Чохральского 164
5.1. Бифуркация осесимметричного стационарного течения в трёхмерное 164
5.1.1. Предварительные замечания и постановка задачи 164
5.1.2. Возникновение стационарных трёхмерных структур 166
5.2. Трёхмерная конвекция в переходном режиме 173
5.2.1. Постановка и вычислительные аспекты задачи 173
5.2.2. Изменение характера течения с ростом числа Релея 175
5.2.3. Влияние вращения ёмкости и центрального тела на конвекцию расплава 187
5.2.4. Осреднённые поля и их сравнение с осесимметричным решением 194
5.2.5. Теплоотдача на интерфейсе кристалл/расплав 198
Заключение 200
Литература 202
- Неустойчивость подогреваемого снизу вращающегося слоя
- Пространственная дискретизация и расчёт невязок
- Развитие неустойчивости центрального вихря: квазипериодический и стохастический режимы
- Влияние вращения ёмкости и центрального тела на конвекцию расплава
Введение к работе
Современный уровень развития полупроводниковой электроники во многом связан с новейшими достижениями в технологии получения электронных материалов. Монокристаллы больших размеров выращивают из расплава, при этом основным методом их получения является метод Чохральского, основу которого составляет процесс кристаллизации полупроводника на охлаждаемой затравке, расположенной на свободной поверхности расплава, помещённого в подогреваемый тигель.
Полученные кристаллы применяются в производстве интегральных схем, размеры которых непрерывно растут, при этом идёт процесс резкого увеличения числа элементов на одной пластине. Из этого проистекает необходимость сокращения числа дефектов и повышения однородности кристаллов при одновременном увеличении их размеров. Гидродинамика и тепломассообмен в жидкой фазе (расплаве) практически полностью определяют свойства получаемых кристаллов, в том числе и в отношении количества дефектов.
Течение расплава при совместном действии ряда факторов, в первую очередь, плавучести в поле силы тяжести из-за перепада температур между поверхностями тигля и кристалла, вращения тигля и вращения кристалла, характеризуется многообразием взаимосвязанных явлений, полное представление о которых можно получить при последовательном моделировании каждого из физических эффектов. Имеющийся на сегодняшний день объём сведений, в том числе фундаментального характера, о свойствах режимов течения полупроводниковых расплавов, характеризующихся малыми значениями числа Прандтля, крайне недостаточен.
Чрезвычайно важным представляется получение детальных знаний о трёхмерной нестационарной конвекции в емкостях, вращающихся вокруг вертикальной оси, в приложении к оптимизации метода выращивания кристаллов полупроводников из расплава (метода Чохральского). Прямое численное моделирование трёхмерной нестационарной конвекции на сегодняшний день является наиболее перспективным инструментом при исследовании движения расплава. Его применение позволит существенно дополнить имеющийся скудный объём экспериментальной информации по существующим устройствам, а в последующем в существенной степени заменить дорогостоящие экспериментальные исследования при отработке новых установок.
Исходя из этих соображений, определены основные цели работы. Диссертационная работа направлена на
1) разработку и реализацию метода численного моделирования трёхмерных квазипериодических и стохастических режимов конвекции в областях сложной геометрии;
2) численное моделирование конвекции в моделях тигля метода Чохральского (рассматриваются цилиндрические и кольцевые ёмкости) и описание проявляющихся в них бароклинной неустойчивости и неустойчивости центрального вихря;
3) прямое численное моделирование конвекции в емкостях с геометрией, типичной для тиглей метода Чохральского, при типичных же тепловых граничных условиях, вплоть до чисел Релея порядка 105, а также на выявление условий бифуркации осесимметричного стационарного течения в трёхмерное стационарное и условий перехода к стохастическому характеру течения.
Первая глава диссертации посвящена обзору работ по экспериментальным исследованиям и численному моделированию конвекции расплава. Отражены основные результаты экспериментального моделирования конвекции как в реальных установках метода Чохральского, так и в моделях тигля упрощённой геометрии, в том числе с использованием модельных жидкостей. На основании имеющихся данных сделан вывод о том, что в условиях, характерных для промышлен-но применяемых в настоящее время устройств выращивания кристаллов по методу Чохральского, течение носит сугубо трёхмерный характер. Доступная к настоящему времени экспериментальная информация как по существующим устройствам, так и по моделям, не позволяет получить полное представление о гидродинамике и теплообмене в расплаве. Проанализированы предпосылки к применению прямого численного моделирования трёхмерной нестационарной конвекции расплава. Сопоставлены имеющиеся на сегодняшний день в литературе результаты расчётов в трёхмерной постановке.
Во второй главе рассматриваются математическая модель и численный метод, на основе которых были проведены расчёты конвекции. Течение расплава описывается системой трёхмерных нестационарных уравнений Навье-Стокса в совокупности с уравнением энергии. Для учёта сил плавучести в поле силы тяжести и центробежной силы используется приближение Буссинеска. При проведении расчётов использовалась новая версия разработанного на кафедре гидроаэродинамики СПбГТУ программного комплекса SINF-2, которая позволяет проводить расчёты стационарных и нестационарных течений несжимаемой жидкости или газа, развивающихся в общем случае в областях сложной геометрии. Используются многоблочные структурированные неравномерные сетки, согласованные с границами области течения. Для получения нестационарных решений реализованы схемы как первого, так и второго порядка. При этом на каждом шаге по физическому времени итерации осуществляются по методу искусственной сжимаемости. Для расчёта конвективных слагаемых применяется схема QUICK, а при дискретизации дифференциальных операторов, отражающих действие вязкости, применяется центральная разностная схема второго порядка.
В третьей главе описаны проведённые методические и тестовые расчёты. С целью отработки методик расчёта нестационарных течений с использованием реализованного метода численного моделирования квазипериодических и стохастических режимов конвекции рассмотрено возникновение автоколебательных режимов в квадратных полостях с разнонагретыми стенками как при умеренных, так и при малых числах Прандтля. Для оценки требуемой размерности сеток при проведении трёхмерных расчётов естественной конвекции выполнено моделирование стационарной ламинарной конвекции в кубической полости с двумя боковыми стенками, поддерживаемыми при различных температурах. Рассмотрено также развитие автоколебательного режима конвекция при малой надкритично-сти в неподвижной цилиндрической модели тигля метода Чохральского с примыкающим неподвижным диском. Для указанных задач выполнено сопоставление с экспериментальными и расчётными данными других авторов.
Четвёртая глава посвящена численному моделированию конвекции в моделях тигля метода Чохральского (рассматриваются кольцевые и цилиндрические ёмкости) и описанию проявляющихся в них бароклинной неустойчивости и неустойчивости центрального вихря. Получены данные о полях скорости и температуры для условий, при которых во вращающихся емкостях возникают структуры, обусловленные бароклинной неустойчивостью; воспроизведены наблюдавшиеся в экспериментах режимы регулярных волн для больших чисел Прандтля и режимы негеострофической турбулентности для малых чисел Прандтля. Выявленные свойства конвекции во вращающихся емкостях и полученные подробные данные по полям скорости и температуры являются важным дополнением к весьма скудной экспериментальной информации по локальным характеристикам и будут способствовать более рациональным постановкам последующих опытов. Также в четвёртой главе поставлена и исследована модельная задача изотермического течения во вращающейся ёмкости с примыкающим противовращающимся центральным диском, направленная на выявление роли неустойчивости центрального вихря в развитии трёхмерных автоколебательных процессов в установках метода Чохральского. Проведённые вычисления показали, что эта роль может быть определяющей и в режимах тепловой конвекции.
В пятой главе рассматривается конвекция во вращающихся емкостях, усложнённая форма которых типична для тиглей систем метода Чохральского. Выявлены условия бифуркации осесимметричного стационарного течения в трёхмерное стационарное течение при умеренных числах Релея в невращающемся тигле. Проведены параметрические расчёты трёхмерных квазипериодических и стохастических режимов конвекции во вращающемся тигле в условиях противо-вращения кристалла. Получен обширный набор данных об эволюции структуры течения при изменении определяющих параметров и спектральном составе пульсаций.
В заключительных выводах представлены основные результаты и рекомендации, полученные в работе.
Неустойчивость подогреваемого снизу вращающегося слоя
Работы по физическому моделированию течения расплава кремния при выращивании монокристаллов по методу Чохральского появились в 50-е годы, что объяснялось бурным развитием полупроводниковой промышленности. Обзор основных этапов экспериментального моделирования гидродинамических процессов в расплаве для метода Чохральского в нашей стране приведён в [4]. В большей части отечественных работ, направленных на экспериментальное моделирование течения расплава и исследование структуры течения, использовались не расплавы полупроводников, а модельные жидкости, например, вода [2] или водоглицериновые смеси [4]. Возможность непосредственного переноса результатов, полученных для модельных жидкостей, на течение расплавов полупроводников, характеризующихся малыми значениями числа Прандтля, требует дальнейших исследований.
И в зарубежных работах первоначально изучались преимущественно закономерности конвекции модельных жидкостей [24], не оставлено это направление и до последнего времени. Параллельно началось проведение экспериментов, в которых непосредственно изучалось течение расплавов полупроводников, в том числе в натурных условиях, особенно бурный рост таких исследований наблюдается с конца 80-х годов. Одним из примеров появившихся в 80-е годы экспериментальных работ, направленных на доказательство существенного проявления трёхмерных эффектов при конвекции расплава в установках метода Чохральского, является статья [91]. В ней описываются эксперименты в лабораторной ростовой установке с тиглем довольно малого размера (Rc = 38 мм, Нс = 37,5 мм, Ra 10 ). Для выявления эффектов плавучести в чистом виде тигель и кристалл были остановлены, вытягивание кристалла на время измерений было прекращено. Наблюдалась непосредственно структура течения путём визуализации траекторий вольфрамовых пробных частиц в рентгеновских лучах. В ходе экспериментов было получено наглядное подтверждение существенной трёхмерности течения, хотя использованная при проведении опытов однолучевая рентгеновская установка, как свидетельствуют более поздние данные, даёт весьма искажённое представление о структуре пространственного течения [155]. Дополнительно исследовался вопрос о степени влияния на трёхмерный характер течения неосесимметричности тепловых граничных условий. Указанный фактор был признан существенным, однако дальнейшего развития в литературе этот вопрос не получил.
В работе Kakimoto et al. [92] описаны результаты более аккуратного эксперимента, проведённого в том же самом тигле, что и в предыдущем случае, но в присутствии вращения (тигель и кристалл вращаются как единое целое). Впервые наблюдался относительно упорядоченный вихрь, практически неподвижный с точки зрения наблюдателя, помещённого во вращающуюся систему координат. В то же время с точки зрения неподвижного наблюдателя течение имело хаотический характер. Как и в предыдущем случае, структура течения наблюдалась путём визуализации траекторий пробных частиц в рентгеновских лучах. Описать поведение трёхмерных структур стало возможным благодаря появлению двухлу-чевой рентгеновской установки (подробное описание установки и методики проведения экспериментов дано в [155]). Экспериментальные данные, дополненные результатами численного моделирования, привели авторов к выводу, что полученное течение является бароклинной волной. Этот вывод был закреплён результатами экспериментов Watanabe et al. [154], проведённых в широком диапазоне определяющих параметров (варьировались перепад температуры и скорости вращения кристалла и тигля). Таким образом, причины возникновения трёхмерных вихревых структур в расплаве, по мнению данной группы авторов, следует отнести на счёт бароклинной неустойчивости.
Одновременно с вышеупомянутыми исследователями предположение о появлении медленно прецессирующих волновых структур, похожих на наблюдавшиеся во вращающихся кольцевых емкостях [68], высказали Kishida et al. [100], которые провели экспериментальные наблюдения за течением расплава в 16-дюймовом тигле с размещённым на свободной поверхности неподвижным диском (измерялись колебания температуры в нескольких точках мониторинга). Было показано, что при увеличении скорости вращения тигля наблюдается некоторое упорядочивание течения, проявляется ведущая частота колебаний, амплитуда пульсаций на периферии тигля возрастает, а величина колебаний в центре тигля остаётся небольшой. Во вращающейся системе координат направление дрейфа бароклинной волны при этом меняется с совпадающего с направлением вращения тигля на противоположное.
Несколько позднее Yi et al. [161] провели измерения температуры расплава во вращающемся цилиндре с размещённым на свободной поверхности диском с целью подтвердить или опровергнуть присутствие бароклинных волн. Ёмкость и диск вращались в одну и ту же сторону с одинаковыми угловыми скоростями. Показано, что с увеличением вращения для рассмотренного значения числа Релея (Ra= 1,8x105) поле температуры становится неосесимметричным. В работе зафиксировано возникновение бароклинных волн при RoT 102, что существенно превышает значения, полученные при экспериментах в кольце [68]. Характерное волновое число m равняется 2.
Крайне интересными с точки зрения проявления механизма бароклинной неустойчивости следует признать проведённые Lee, Chun [111] эксперименты по измерению колебаний температуры в наполненном ртутью цилиндре с нагретыми боковыми стенками, на свободной поверхности размещён неподвижный диск. Условия эксперимента полностью описаны, что позволяет использовать результаты работы в качестве теста при численном моделировании. Существенным положительным моментом исследования следует признать также тот факт, что вместо горячего расплава кремния использовалась ртуть. В отличие от опытов с расплавом кремния, где всегда есть существенная неопределённость при фиксации характерного перепада температуры в силу невозможности прямого контроля за температурой стенок тигля, при проведении экспериментов в ртути возможен точный контроль тепловых граничных условий. В [111] для разных геометрических параметров (А= 1...2,5, (3=1) продемонстрировано появление регулярных волновых структур с волновыми числами m от 1 до 3, обладающих характерными свойствами бароклинных волн (как то: скорость прецессии волны, однородность волны в вертикальном направлении и т.д.). В то же время, были зафиксированы и существенные отличия в поведении волновых структур в цилиндрической ёмкости по сравнению с вращающимся кольцом, в частности, значения волновых чисел в области регулярных волн для цилиндра существенно меньше значений для кольца, где максимальное значение m как раз наблюдается в области наибольшей регуляризации течения. Это расхождение, отмеченное и в [100], по утверждению авторов следует отнести на счёт различий в геометрии рассмотренных емкостей.
Пространственная дискретизация и расчёт невязок
Показано, что учёт вращения диска и ёмкости очень сильно меняет структуру течения. Дополнительно в указанной работе исследовалось влияние на течение расплава термокапиллярных эффектов. Показано, что для данных условий эффект Марангони не слишком сильно влияет на структуру течения и поле температуры, но для получения количественных результатов его желательно учитывать. Подробные сведения о нескольких рассмотренных режимах течения (варьировались высота ёмкости и скорости вращения ёмкости и диска) приведены в [151].
Все остальные результаты численного моделирования во многом носят качественный характер, позволяющий лишь в какой-то степени прояснить структуру течения и выявить основные механизмы неустойчивости. Основное внимание при проведении исследований уделялось проявлениям бароклинной неустойчивости. Так, направленные на исследование влияния тепловых граничных условий и вращения на поведение расплава численные эксперименты Kakimoto et al. [92, 93] продемонстрировали возникновение вихревых структур, природа которых была отнесена к бароклинной неустойчивости. Отличительной особенностью описываемых расчётов явился тот факт, что в качестве граничных условий на стенках ёмкости были взяты распределения температуры, имитирующие полученные в ходе решения осесимметричной задачи о глобальном теплообмене в установке метода Чохральского для тигля реальной геометрии. Расчёты проводились с использованием коммерческого кода FLUENT. Получено качественное совпадение с экспериментальными результатами, опубликованными в [92]. Yi et al. [161] провели численное моделирование течения расплава во вращающемся цилиндре, на свободной поверхности которого размещён диск, вращающийся в ту же сторону и с той же скоростью, что и ёмкость. Расчёты проводились при значениях числа Релея, находящихся в диапазоне 103-105 (Ro = 0,7...23). Боковые стенки и дно поддерживались при одной и той же постоянной температуре, что привело к обычному на практике сосуществованию вертикального и горизонтального градиентов температуры. При сравнительно медленном вращении в работе рассмотрены режимы, в которых проявляется конвекция, свойственная подогреваемому снизу вращающемуся слою, произведена оценка критического числа Релея для разных скоростей вращения и проведено сопоставление с аналитическим решением [8, 59] для безграничного слоя, полученным на основе линейной теории устойчивости. Критические значения числа Релея для рассмотренного Yi et al. случая оказались значительно выше аналитической кривой практически во всём диапазоне угловых скоростей вращения. Для сравнительно быстрого вращения ёмкости впервые показано образование системы вихрей на периферии расплава, прецессирующей в направлении вращения ёмкости. На основе аналогии с результатами [68] для конвекции ртути во вращающемся кольце был сделан вывод о существенном проявлении бароклинной неустойчивости при быстром вращении.
В работе [138] численное моделирование течения расплава во вращающейся цилиндрической ёмкости без диска было проведено для Ra= 10 при варьировании числа Россби. Сопоставление результатов численного моделирования с приведёнными в этой же работе результатами экспериментов для Ra«2xl06, показало качественное совпадение. В работе приведены распределения скорости и температуры, подтверждающие экспериментальные данные и позволяющие наглядно судить о появлении прецессирующих волновых структур с волновыми числами 3-5. Как и в эксперименте, был выявлен переход от прецессии "с опережением" при слабом вращении к прецессии "с отставанием" при сильном вращении.
Наиболее приближенными к реально достигаемым значениям числа Релея стали расчёты Tanaka et al. [143], также проведённые в цилиндрической ёмкости в отсутствие диска. Были представлены результаты расчётов для трёх скоростей вращения (Ro = 0,5 - 16), которые показали качественное согласование с опубликованными в этой же работе распределениями температуры на свободной поверхности, полученными с помощью CCD-камеры. Авторы не стали относить причины проявления неустойчивости на счёт бароклинной неустойчивости.
Расчёт конвекции при таких же условиях для одной скорости вращения (Ro = 0,5) проведён Basu et al. [49]. В работе обсуждается проблема выбора сеток, на конкретных примерах показаны преимущества многоблочного подхода. Решение в определённой степени согласуется с результатами [143]. Показано, что при сильном вращении течение имеет ячеистую структуру, связанную с присутствием конвективных столбиков, одновременно на периферии ёмкости присутствуют слабые волны, прецессирующие в направлении вращения ёмкости. На основании приведённых результатов очевидно, что даже использованной мелкой сетки размерностью 337 тыс. ячеек недостаточно для количественного описания течения при столь высоком значении числа Релея без применения подсеточных моделей.
Все указанные выше работы были посвящены численному моделированию конвекции жидкости с малым значением числа Прандтля, характерным для расплава кремния. Однако область применения метода Чохральского достаточно широка и не ограничивается лишь выращиванием монокристаллов кремния. В частности, по методу Чохральского выращивают так называемые оксидные кристаллы (LiNb03, AI2O3 и т.д.), расплавы которых характеризуются сравнительно большими значениями числа Прандтля порядка 10. Как показали результаты экспериментального и численного моделирования [86], движение расплава при этом также является автоколебательным и существенно трёхмерным, что, однако, вызвано принципиально другими механизмами по сравнению с конвекцией расплава кремния, прежде всего, в силу отсутствия вращения тигля. Течение в очень сильной степени определяется влиянием термокапиллярной конвекции, в результате которого на поверхности расплава возникает так называемая "лепестковая", или "спицевидная" структура ("spoke pattern").
Развитие неустойчивости центрального вихря: квазипериодический и стохастический режимы
При проведении расчётов использовалась новая версия разработанного на кафедре гидроаэродинамики СПбГТУ программного комплекса SINF-2; детальное описание алгоритма исходной версии программы дано в [141]. Программный комплекс позволяет проводить расчёты стационарных и нестационарных течений несжимаемой жидкости или газа, развивающихся в общем случае в областях сложной геометрии. Используются многоблочные структурированные неравномерные сетки, согласованные с границами области течения. Уравнения движения записаны относительно декартовых компонент скорости.
К настоящему времени разработано значительное число достаточно совершенных методов численного решения уравнений Навье-Стокса. Известно, что при численном решении уравнений Навье-Стокса в приближении Буссинеска наибольшие трудности возникают из-за слабой связи между уравнениями движения и неразрывности. Уравнения движения, в частности, выражают связь между полем давления и полем скорости, подстраивая последнее к получающемуся полю давления. В то же время для случая постоянной плотности уравнение неразрывности (баланса массы) не содержит члена, зависящего непосредственно от поля давления. Эффективным способом введения исходно отсутствующей связи между полями скорости и давления в уравнение баланса массы является добавление к нему члена, моделирующего искусственную сжимаемость и содержащего производную от давления по фиктивному времени, не связанному с реальным, физическим временем процесса. В случае стационарного течения модифицированное уравнение баланса массы не имеет физического смысла до тех пор, пока в процессе вычисления не будет достигнуто стационарное состояние потока.
Метод искусственной сжимаемости был почти одновременно независимым образом предложен рядом исследователей [5, 61]. С тех пор он получил достаточное развитие и с успехом использовался во многих работах для решения задач несжимаемых стационарных течений. Пейре и Тейлор [29] представили детальное обсуждение метода и обзор численных схем. Решение уравнений Навье-Стокса для несжимаемой жидкости в программном комплексе SINF-2 осуществляется по методу искусственной сжимаемости.
Дискретизация пространственных операторов дифференциальных уравнений сохранения в программном комплексе выполнена по методу контрольного объёма со вторым порядком точности. В [41] дается обзор этого метода, вводится он как частный случай метода взвешенных невязок. Основная идея метода состоит в следующем [27]: расчётную область разбивают на непересекающиеся целиком заполняющие область контрольные объёмы таким образом, что каждая узловая точка содержится в одном контрольном объёме. Дифференциальное уравнение, отвечающее некоторому закону сохранения, интегрируют по каждому контрольному объёму, используя теорему Гаусса. Для вычисления интегралов используют кусочные профили, которые описывают изменение подынтегральной функции между узловыми точками.
Дифференциальное уравнение выражает закон сохранения для бесконечно малого объёма, а полученный в результате интегрирования его дискретный аналог точно выражает закон сохранения для конечного контрольного объёма. Таким образом, решение удовлетворяет точному интегральному балансу при разбиении расчётной области на любое количество объёмов, а значит и при использовании грубых сеток, то есть метод обладает хорошими консервативными свойствами
Значения искомых величин определяются в центрах расчетных ячеек (контрольных объёмов). Для расчёта конвективных слагаемых применяется противо-поточная разностная схема высокого порядка, обобщающая схему QUICK и линейную противопоточную схему второго порядка [115, 140]. Для дискретизации дифференциальных операторов, отражающих действие вязкости, применяется центральная разностная схема второго порядка. Расчеты осуществляются на совмещённой сетке, что может привести к нефизическим осцилляциям давления, для их подавления используется специальная интерполяционная процедура коррекции для вычисления массовых потоков в уравнении неразрывности [133].
Продвижение по фиктивному времени осуществляется с помощью неявного метода первого порядка, основанного на линеаризации определяемого исходными уравнениями дифференциального оператора и на методе приближенной факторизации [50]. Система алгебраических уравнений, полученная в результате использования неявного метода приближенной факторизации, решается методом матричной прогонки для нахождения компонент вектора скорости и скалярной прогонки для остальных переменных.
В используемой для расчётов многоблочной версии программного комплекса SINF-2 на каждом шаге по фиктивному времени осуществляется обмен значениями рассчитываемых переменных вдоль поверхностей, по которым происходит стыковка блоков. При этом используется концепция вспомогательного виртуального блока, что обеспечивает полную "прозрачность" межблочных границ и сохранение консервативных свойств разностной схемы.
Заложенные в комплекс возможности были верифицированы в ходе численного моделирования целого ряда задач гидродинамики и теплообмена: моделирование течений в вихревых устройствах со сложной геометрией [16, 38, 43], вращающихся кольцевых полостях [17, 96] и каналах [37, 94, 97, 98, 142]. С помощью программ SINF и SINF-2 решались также некоторые задачи промышленной аэродинамики [6, 15, 19, 35, 46, 48].
В большинстве случаев расчётная область при решении уравнений Навье-Стокса имеет сложную форму, ее границы не совпадают в физическом пространстве с координатными линиями ортогональных систем координат (например, декартовой или цилиндрической), что требует применения интерполяции на сеточные линии при реализации граничных условий. Это ведет к усложнению расчётного алгоритма и локальной потере точности численного решения в приграничных областях. Указанные трудности преодолеваются введением согласованных с границами расчётных областей криволинейных координатных систем [42] .
Влияние вращения ёмкости и центрального тела на конвекцию расплава
В качестве исходного теста для верификации программного кода и отработки методик численного моделирования стационарных и нестационарных термоконвективных течений были выбраны два частных случая ставшей уже канонической задачи о естественной конвекции в прямоугольной полости с разнона-гретыми вертикальными стенками [36]. В литературе представлена подробная экспериментальная и расчётная информация о свойствах течения в широком диапазоне определяющих параметров, к которым относятся число Релея Ra, число Прандтля Рг и коэффициент формы А. В продолжение последних двух десятилетий рассматриваемая задача, имеющая широкое прикладное значение, стала одним из наиболее популярных тестов для программных кодов, используемых при расчёте термоконвективных течений. В настоящей работе рассматривались стационарные и нестационарные течения в квадратной полости (А= 1) при Рг = 0,71 (воздух) и Рг = 0,023 (расплав галлия).
Обзор посвященных изучению стационарной конвекции воздуха в квадратной полости работ, дающий достаточно полную картину состояния проблемы на начало 90-х годов, содержится в [47]. Основа постоянно пополняющейся базе численных решений данной задачи была положена De Vahl Davis в работе [62], где исследовалось стационарное ламинарное течение воздуха в квадратной каверне при числах Релея в диапазоне 103 - 106. Автором был предложен набор образцовых решений, погрешность которых не превышала 1%. Стационарные ламинарные решения в том же диапазоне чисел Релея были получены Markatos, Pericleous [118]. Кроме того, в указанной работе для расчёта турбулентных режимов конвекции в квадратной полости для чисел Релея вплоть до 10 впервые применялась k-є модель турбулентности.
В последние годы усилия исследователей были направлены на моделирование переходных режимов течения в квадратной полости и изучение механизмов неустойчивости. Основа этому была заложена в начале 80-х годов Patterson, Imberger в работе [128], где впервые было предсказано существование периодического режима течения. Различные подходы к изучению периодических и квазипериодических режимов течения описаны Janssen, Henkes в работе [84], там же приведены основные результаты моделирования перехода в квадратной полости, полученные в 80-х - первой половине 90-х годов, с соответствующими ссылками. Было установлено [84], что в случае конвекции воздуха в квадратной полости в ходе перехода от ламинарного течения к турбулентному реализуется сценарий Рюэля и Такенса [136]. При достижении некоторого критического числа Релея порядка 108 стационарное течение становится периодическим. Затем, при дальнейшем росте числа Релея, происходит бифуркация к квазипериодическому режиму с двумя ведущими частотами, который постепенно вырождается в хаотический режим. Расчёты Paolucci, Chenoweth [127] продемонстрировали, что первая бифуркация связана с потерей устойчивости течения в угловых областях при повороте воздуха, движущегося в вертикальных пограничных слоях. Эта неустойчивость, как показано в [132], носит термоконвективный характер. В [127] показано также проявление второй бифуркации, в ходе которой проявляется неустойчивость вертикальных пограничных слоев вдоль изотермических стенок. В работе [84] расчёты проводились в диапазоне чисел Релея вплоть до 7,5х108, причём при Ra = 2x108 был получен периодический режим, при Ra = 3x108 - квазипериодический режим, а при Ra = 7,5xl08 - хаотический режим. Подробная экспериментальная информация о течении воздуха в квадратной полости при хаотическом режиме, Ra = 1,58x109, содержится в работах Tian, Karayiannis [144, 145], где приведены как осреднённые, так и пульсационные характеристики течения.
Особенностью конвекции в прямоугольных областях при малых числах Прандтля является тот факт, что уже при относительно небольших значениях числа Релея наблюдается переход к нестационарным режимам течения. Поэтому уже с конца 70-х годов при лабораторном и численном моделировании стационарные режимы сами по себе практически не рассматривались [45], в основном изучались периодические и квазипериодические режимы течения, а также осуществлялся поиск критических значений параметров, определяющих появление колебаний (см., например, [7]). Специфика изучения конвекции при малых числах Прандтля заключается ещё и в том, что основное внимание исследователей привлекали мелкие полости. Это объясняется тем, что конвекция в мелких полостях с разнонагретыми вертикальными стенками моделирует процесс выращивания кристаллов из расплава по горизонтальному варианту метода Бриджмена и метода плавающей зоны. Так, в классической экспериментальной работе Hurle et al. [83], показавшей существенное влияние трёхмерных эффектов при нестационарном режиме течения расплава галлия, исследовалось возникновение осцилляции именно в мелких полостях, рассматривались полости со значением А = 2 и выше. На симпозиуме GAMM [130] была подробно проанализирована задача конвекции при А = 4. Данная постановка была провозглашена основой для построения базы данных образцовых решений ("benchmark solution") для тестирования алгоритмов при решении задач конвекции. За прошедшее десятилетие этот частный случай конвекции подвергся детальному анализу множеством авторов [10, 45, 51, 52, 71, 72,119,124,125].
Численному моделированию конвекции жидкости в квадратной полости при малом значении числа Прандтля посвящено сравнительно небольшое число работ. Гельфгат и Мартузан [7] провели оценки критических значений числа Ре-лея и рассчитали колебательные режимы при Pr = 0.02, расчёты носили лишь качественный характер в силу того, что при применении метода Галёркина использовалось малое число базисных функций. Течение в квадратной каверне при очень малых значениях числа Прандтля (Рг = 0,001, Рг = 0,005 и Рг = 0,01) рассматривалось в работе Mohamad, Viskanta [122], где приведены сравнительно грубые оценки критических значений числа Релея. Подробное исследование перехода к периодическим режимам течения проделано Bergman, Ball в [53]. Проведённый объём вычислений позволил определить критические значения числа Ре-лея при 0,005 Рг 0,035. В работе применялся метод Чебышёва. Gelfgat et al. [71] получили диаграмму устойчивости течения в полости со свободной поверхностью в диапазоне 1 А 10 при Рг = 0 и Pr = 0,015. В работе [160], посвященной исследованию пригодности схемы QUICK для решения ряда задач гидродинамики и теплообмена, задача о нестационарной конвекции в квадратной полости при Рг = 0,025, Ra = 5х105 рассматривалась в качестве одного из тестовых примеров. В работе приведена картина эволюции среднего числа Нуссельта на изотермической стенке.
