Содержание к диссертации
Введение
Глава 1 . Использование экспоненциальных схем в алгоритмах решения уравнений Навье-Стокса и конвективного теплопереноса ., 18
1.1. Экспоненциальная схема для решения уравнений Навье-Стокса и конвективного теплопереноса 21
1.2. Результаты расчетов двумерных тестовых задач 27
1.3. Использование экспоненциальных схем для решения трехмерных задач, описываемых уравнениями Навье-Стокса 37
1.4. Выводы 48
Глава 2. Исследование конвекции раствор-расплава при выращивании кристаллов в случае неоднородного разогрева боковых стенок ростового тигля 49
2.1. Математическая модель ростовой установки для получения кристаллов методом Чохральского при неоднородном разогреве боковых стенок тигля 51
2.2. Алгоритм численного решения 54
2.3. Анализ температурного поля и структуры течений в условиях циклически изменяющегося разогрева боковых стенок ростового тигля 65
2.4. Анализ условий роста кристалла при стационарном неоднородном разогреве боковых стенок ростового тигля 69
2.5. Выводы 74
Глава 3. Процессы тепло- и массонереноса при гидротермальном росте кристаллов 76
3.1. Математическая модель процессов конвективного тепло- и массонереноса в автоклаве 78
3.2. Описание численного алгоритма 81
3.3. Численные исследования процессов гидротермального синтеза при получении искусственных изумрудов 86
3.4. Выводы ...98
Заключение 99
Список литературы
- Результаты расчетов двумерных тестовых задач
- Использование экспоненциальных схем для решения трехмерных задач, описываемых уравнениями Навье-Стокса
- Анализ температурного поля и структуры течений в условиях циклически изменяющегося разогрева боковых стенок ростового тигля
- Численные исследования процессов гидротермального синтеза при получении искусственных изумрудов
Введение к работе
Настоящая работа посвящена численному моделированию и исследованию структур течений, процессов тепло- и массопереноса в установках по выращиванию монокристаллов, а также развитию и разработке алгоритмов реализации возникающих при этом задач.
Получение монокристаллов является сложной проблемой, так как на процесс влияет множество различных факторов, связанных с конвекцией, физическими и химическими свойствами веществ, особенностями технологических установок [21]. Вследствие повышения требований к качеству кристаллов для современной электроники и лазерных установок, становится важным управление ростовыми процессами, напрямую зависящими от конвективного тепломассопереноса.
Принципиальные трудности в управлении процессами гидродинамики, тепло- и массообмена при росте кристаллов вызваны тем, что все они взаимосвязаны. В то же время на практике во многих случаях достаточно иметь представление о структуре существующих течений, распределении температурных полей и их временных характеристиках. Изучение процессов в такой постановке является актуальным при разработке новых технологий ростовых процессов.
В настоящее время самое широкое распространение получили
методы математического моделирования изучаемых процессов, так как
проведение физических экспериментов с использованием реальных
ростовых установок требует больших материальных затрат и
сопровождается технологическими трудностями. Сложность
рассматриваемых задач требует создания двух- и трехмерных моделей с использованием уравнений Навье-Стокса. В этих условиях для получения решений с нужной точностью используется аппарат численных методов с применением ЭВМ.
На сегодняшний день численное моделирование с использованием уравнений Навье-Стокса сформировалось как самостоятельное направление в механике жидкости, ее приложениях к гидродинамике, машиностроению, энергетике, а также к изучению природных явлений. Постоянные требования все более точных расчетов характеристик рабочих процессов при поиске оптимальных конструкторских и технологических решений делают необходимым развитие математических моделей, а также алгоритмов их реализации [22, 24, 78].
Целью диссертационной работы является исследование возможности использования экспоненциальных разностных схем в алгоритмах численной реализации уравнений Навье-Стокса и конвективного теиломассопереноса, а так же применение таких алгоритмов для анализа структуры течения и температурного поля в жидкости при получении кристаллов в ростовых установках.
Новизна исследований состоит в следующем;
Проведен анализ применения экспоненциальных разностных схем, полученных после преобразования уравнений Навье-Стокса к самосопряженному виду, в алгоритмах численной реализации многомерных задач.
Исследована структура пространственного (трехмерного) течения и процессов конвективного тепломассопереноса в раствор-расплаве при получении кристаллов методом Чохральского в условиях стационарного и нестационарного азимутально-неоднор одного нагрева боковых стенок ростового тигля.
Впервые выполнено комплексное численное исследование пространственных полей температуры, структуры течения и концентрации растворенного вещества в установке гидротермального роста кристаллов при описании системы в виде цилиндрического сосуда, заполненного жидкостью и частично пористым слоем, при наличии затравки.
В первой главе диссертационной работы проводится анализ использования экспоненциальных схем в алгоритмах численной реализации уравнений Навье-Стокса и конвективного тепломассопереноса.
При описании поведения вязкой жидкости используются уравнения Навье-Стокса. Впервые эти уравнения были сформулированы Навье в 1821 году как обобщение уравнений Эйлера, полученных в 1751 году. Позднее, в 1845 году, Стоксом была представлена их окончательная запись. Однако до сих пор отсутствует общий аналитический метод получения решений уравнений Навье-Стокса.
В моделях, рассматриваемых в работе, предполагается, что жидкость вязкая и несжимаемая, что ее плотность зависит только от температуры и не зависит от давления, другие физические свойства вещества постоянны для рассматриваемого объема жидкости. Эти условия наряду с предположением о малости отклонений системы от состояния гидростатического равновесия составляют основу приближения Обербека-Буссинеска [5,12].
На сегодняшний день имеется большое количество методов решения уравнений Навье-Стокса, применяемых, главным образом, для решения двумерных задач, Чаще всего при этом используются переменные вихрь скорости - функция тока [3,18,22,24,38,40,64,85]. Однако в случае расчетов пространственных течений их описание через вихрь скорости и функцию тока приводит к большому числу зависимых переменных. К тому же постановка граничных условий для вихря скорости на твердой поверхности затруднительна. Поэтому, как правило, расчеты пространственных течений несжимаемой жидкости проводятся в простейших переменных - давление, составляющие скорости [11,47,56,58]. Существующая при этом проблема определения давления может быть разрешена одним из двух основных способов.
Во-первых, давление может вычисляться из уравнения Пуассона с
использованием прямых или итерационных методов отдельно от
уравнений движения [24]. Другой способ заключается в том, что давление рассчитывается одновременно со скоростью. Например, метод искусственной сжимаемости, впервые описанный в [46], суть которого заключается в добавлении производной по времени от давления в уравнение неразрывности. При этом невязкая часть уравнений Навье-Стокса становятся гиперболической по времени и для нахождения стационарных решений модифицированной системы можно применять метод установления. Распространение данного подхода на нестационарный случай возможно с помощью так называемого "двойного шага по времени" [43,61]. Наряду с этим, для определения давления одновременно со скоростью в последние годы все более широкое применение находят так называемый алгоритм Узавы [41,44,45,53] и метод Gauge [16,48].
Исторически после появления вычислительных машин методы конечных разностей (МКР) были наиболее распространенным аппаратом решения задач математической физики. Однако впоследствии пальма первенства в теоретических исследованиях и практическом применении (особенно за рубежом) перешла к методам конечных элементов (МКЭ), основанном на аппроксимации вариационных постановок, что отражается в количестве статей и монографий (см. обширную библиографию в [1]).
За последние десятилетия активизировались исследования по методам конечных объемов (МКО) [65], или балансным методам. В определенном смысле эти алгоритмы сближают технологии методов конечных разностей и конечных элементов. Методологически же они представляют собой «новую волну» балансных, или консервативных, разностных схем, основанных на приближениях интегральных соотношений, являющихся следствием дифференциальных уравнений. Метод баланса (интегро-интерполяционный метод) был предложен Самарским А.А. в работе [28] в конце 50-х годов и активно используется в вычислительной практике [9,17,25].
Применению методов конечных разностей для решения задач механики жидкости и газа посвящена обширная литература. Упомянем наиболее известные монографии российских и зарубежных авторов С.К. Годунова и B.C. Рябенького [6], Р.Д. Рихтмайера и К.Мортона [23], В. Вазова и Дж. Форсайта [2], Г.И. Марчука [13-15] и А.А. Самарского [26,27,29-36] с соавторами, К. Флетчера [38], П. Роуча [24].
Как правило, сеточные методы обосновываются и применяются к численному решению дифференциальных задач, относительно которых известны существование, единственность, корректность и гладкость решения. Однако иногда вычисления приходится проводить в таких сложных случаях, когда теоретические вопросы являются открытыми. Тогда адекватность расчетов может основываться на методически грамотно построенных вычислительных экспериментах, когда разностные решения на последовательности сгущающихся сеток подтверждают практически и устойчивость, и сходимость результатов к предельным значениям.
На сегодняшний день нашли применение два направления при построении разностных уравнений. Одно из них заключается в раздельной аппроксимации первой и второй производных. При этом первая производная аппроксимируется направленными разностями против потока или центральными разностями, а вторая производная аппроксимируется стандартным оператором на трехточечном шаблоне [6,19,35]. Другое направление заключается в совместной аппроксимации первой и второй производных. Здесь схема строится после приведения уравнения к дивергентному виду, когда вводится суммарный поток, складывающийся из конвективного и диффузионного потоков [9], либо получение схемы обусловлено преобразованием исходного уравнения к самосопряженному виду[20,59].
Известно, что замена первых производных центральными
разностями позволяет получить аппроксимацию второго порядка. Однако
такой подход непригоден при аппроксимации на трехточечном шаблоне второй производной с малым параметром при ней, так как для получения монотонной схемы существуют ограничения на шаг пространственной сетки [8]. Избежать трудностей позволяет использование направленных против потока разностей для аппроксимации первой производной. Однако, односторонние разности обеспечивают на равномерной сетке аппроксимацию лишь первого порядка. При этом появляется так называемая схемная вязкость, которая может превысить вязкость, описываемую дифференциальным уравнением [8].
Многие исследования посвящены разработке разностных схем, соединяющих достоинства схем с центральными разностями и схем с направленными разностями - второй порядок аппроксимации плюс безусловная монотонность [31]. Основным методом построения монотонных разностных схем для задач, описываемых уравнениями Навье-Стокса и конвективного тепло- и массопереноса, стала их регуляризация. Принцип регуляризации заключается в том, что сначала для исходного уравнения строится простейшая разностная схема, не обладающая необходимыми свойствами аппроксимации, монотонности, устойчивости и т.д. Улучшение качества разностной схемы в необходимую сторону достигается за счет возмущения (регуляризации) коэффициентов, после чего она приводится к некой канонической форме, свойства которой исследованы [27].
Один из способов регуляризации, когда в качестве возмущающих функций используются экспоненциальные функции, называется методом экспоненциальной подгонки [7]. Первые исследования и применение полученных таким образом схем, имеющих коэффициенты в виде экспоненциальных функций, связано с именами A.M. Ильина [8,62], Шарфеттера и Гуммеля [63].
В настоящей работе рассматриваются схемы, получение которых
обусловлено преобразованием исходных уравнений (записанных как в
дивергентной форме, так и в недивергентной) к самосопряженному виду [20,59]. Получаемые после аппроксимации преобразованных уравнений разностные схемы с учетом вида их коэффициентов называются экспоненциальными [31]. Достаточно подробно теоретический анализ таких аппроксимаций приведен в работе А.А. Самарского и П.Н. Вабищевича [31]. К их достоинствам надо отнести то, что на равномерной пространственной сетке они аппроксимируют уравнения Навье-Стокса в приближении Буссинеска и конвективного теплопереноса со вторым порядком точности, а в случае использования неявного алгоритма являются монотонными при любом соотношении пространственных и временных шагов. Однако широкого применения на практике схемы такого вида пока не получили, и их использование требует своего дальнейшего изучения и анализа.
В первом параграфе настоящей главы структура предлагаемого алгоритма представлена на примере одномерной модельной задачи. После преобразования уравнений Навье-Стокса и конвективного теплопереноса, записанных в недивергентном виде, к самосопряженному виду и дальнейшей их аппроксимации можно получить разностную схему, которая является монотонной независимо от того, положительное или отрицательное значение имеет скорость, а также при любом соотношении пространственных и временных шагов. Аппроксимация уравнений Навье-Стокса записанных в дивергентном виде и преобразованных к самосопряженному виду, позволяет получить систему с диагональным преобладанием по столбцам. Показано, что на равномерной пространственной сетке рассмотренные разностные схемы имеют порядок
точности 0(x,h ). Проведено сравнение численного решения с точным на
примере одномерной задачи на основе уравнения Бюргерса.
Предложенный алгоритм был апробирован при решении ряда тестовых задач. Во втором параграфе представлен результат, полученный
при реализации двумерной тестовой задачи о конвекции жидкости со свободной поверхностью в полости прямоугольного сечения. Рассматривается решение широко известной тестовой задачи о конвективном переносе тепла в полости прямоугольного сечения при подогреве сбоку [22,37], полученное с помощью предложенного алгоритма. Дискретизация системы уравнений проводилась на разнесенных равномерных пространственных сетках (по аналогии с методами типа MAC и SIMPLE [52,114]). Для реализации системы линейных алгебраических уравнений применялись итерационные методы.
Численное решение уравнений конвективного теплопереноса является составной частью реализации моделей с использованием уравнений Навье-Стокса в приближении Буссинеска. Поэтому во втором параграфе демонстрируются результаты решения еще одной модельной задачи описываемой двумерным уравнением конвективного теплопереноса в полярных координатах. Эта же задача была решена с использованием треугольных сеток. Использование экспоненциальных схем при численной реализации задач на треугольных сетках упрощает аппроксимацию уравнений и дает возможность получать результаты с высокой точностью.
В третьем параграфе исследуется использование экспоненциальных схем для решения трехмерных задач, описываемых уравнениями Навье-Стокса. Рассматривается модельная задача о конвекции жидкости со свободной поверхностью в цилиндрическом сосуде.
Проводится анализ применимости экспоненциальных схем при
реализации моделей, использующих трехмерные уравнения. С этой целью
рассматривается реализация описанного алгоритма в случае решения
тестовой задачи о конвекции жидкости со свободной поверхностью в
объеме, имеющем форму куба. Проведенные вычислительные
эксперименты подтверждают возможность использования
экспоненциальных схем для численной реализации уравнений Навье-Стокса в приближении Буссинеска, описывающих ламинарные течения.
Вторая глава посвящена исследованию новых технологий управления процессами тепло- и массопереноса при получении монокристалла методом Чохральского. Изделия и элементы, изготовленные из монокристаллов, применяются в качестве различных преобразователей в радиоэлектронике, квантовой электронике, акустике, вычислительной технике и др. Первоначально в технике использовались природные монокристаллы, однако их запасы ограничены, а качество не всегда достаточно высоко. В то же время многие ценные свойства были найдены только у синтетических кристаллов. Поэтому появилась необходимость искусственного выращивания монокристаллов.
Поиску способов управления конвективным тепломассопереносом в ростовых установках во всем мире уделяется большое внимание, поскольку именно эти процессы во многом определяют возможность получения кристаллов с заданными свойствами. Все способы можно разделить на две категории: контактные и бесконтактные. Первые реализуются путем механического вращения кристалла и/или тигля, установкой различных перегородок, мешалок, формообразователен и т. п. [71,76,78,88,89]. К контактным способам можно также отнести воздействие низкочастотных (десятки герц) вибраций, в результате которых в объеме жидкости возникают макротечения, что позволяет подавлять температурные осцилляции на фронте кристаллизации и контролировать распределение примесей в кристалле [67,87].
Вторая категория - бесконтактные способы - основывается на
воздействии физических полей на процессы тепломассопереноса.
Таковыми могут быть: гравитационное поле, магнитное
(электромагнитное) и тепловое. Исследование влияния гравитационного
поля на образование кристаллов связано, главным образом, с космическим
материаловедением, которое выделилось в самостоятельную область
научных исследований. Вместе с тем эксперименты показали, что
существующие на борту космических аппаратов нерегулярные
микроускорения, связанные с работой двигателей ориентации и стабилизации, а также с жизнедеятельностью экипажей, существенно влияют на процессы кристаллообразования, приводя к неоднородности получаемых изделий [69,70,94]. Активно ведутся работы по экспериментальному исследованию и моделированию воздействия магнитных (электромагнитных) полей на конвективные процессы при выращивании кристаллов. Однако такой подход справедлив лишь в случае электропроводящих сред [68,75].
Несколько лет назад А.Е. Кохом была предложена, численно и экспериментально обоснована оригинальная идея выращивания кристаллов в азимутально распределенных стационарных и вращающихся (Heat Field Rotation Method - HFRM) тепловых полях. Суть HFRM заключается в разогреве ростового тигля вертикальными нагревательными элементами, равномерно расположенными по периметру его боковой стенки [72,73,90,91]. Последовательное переключение тепловых источников изменяет температурное поле в жидкости и способствует возникновению в ней азимутальных течений, усиливая перемешивание. Действенность такого подхода демонстрируют результаты, полученные при выращивании кристаллов CsLiB6Oio в циклически меняющемся тепловом поле [90]. Отказ от механического вращения кристалла открывает возможность частичной или полной герметизации ростового пространства, повышает стерильность процесса и избавляет установку от возможных вибраций.
В работе [80] рассмотрена структура течений в цилиндрическом сосуде при циклически меняющемся неоднородном разогреве его боковых стенок в приложении к технологии выращивания кристаллов. Результатами расчетов подтверждена возможность применения бесконтактного метода для гомогенизации раствор-расплава, что значительно упрощает технологию, так как отпадает необходимость переоснащения ростовой установки в разогретом состоянии - замена мешалки на затравкодержатель.
В [74,92,93] приводятся результаты численного моделирования нестационарной термогравитационной конвекции при получении монокристалла методом Чохральского в условиях циклически меняющегося температурного поля в расплаве. На основе анализа данных, полученных в ходе расчетов, сделаны выводы о существовании режимов теплового поля, обеспечивающих максимальное перемешивание кристаллизационной среды при Рг«1. Влияние HFRM на процессы тепломассопереноса в расплавах может способствовать разрушению диффузионного слоя у поверхности кристалла и уменьшению неоднородности распределения примеси в получаемых изделиях.
Вместе с тем особенности конвекции и теплопереноса при использовании HFRM в процессе выращивания монокристалла требуют детального изучения. Проводить достаточно подробные исследования взаимосвязанных процессов конвективного тепломассообмена в расплавах в условиях реальной технологии сложно или невозможно. Это обусловливает, так как течения пространственные, необходимость разработки соответствующей трехмерной математической модели, алгоритма ее решения, позволяющего получать решения, хорошо совпадающие с экспериментальными данными, и проведение вычислительных экспериментов в широком диапазоне меняющихся параметров [79,83].
В первом параграфе второй главы рассматривается математическая модель применения HFRM при получении монокристалла методом Чохральского. Во втором параграфе описан алгоритм ее реализации с использованием экспоненциальной схемы. Численное моделирование процессов конвекции во время роста кристалла основано на совместном решении нестационарных трехмерных уравнений Навье-Стокса б приближении Буссинеска и теплопереноса при заранее определенных условиях нагрева на границе расчетной области.
В третьем и четвертом параграфах приведены результаты, полученные при исследовании зависимости конвективных течений в кристаллизационной среде, обладающей свойствами Bii2Si02o (BSO) с высокой вязкостью (Рг-26), от параметров внешнего теплового поля [77,81].
Третья глава диссертации посвящена исследованию процессов тепло- и массопереноса при гидротермальном росте кристаллов. Гидротермальный рост широко применяется при промышленном выращивании кристаллов [107,111,112]. В этом процессе используется водный раствор под высоким давлением при температуре, необходимой для растворения и рекристаллизации веществ, нерастворимых при обычных условиях. Привлекательность такой технологии заключается в том, что температура, при которой проходит рабочий процесс, сравнительно не высока, и поэтому растущий кристалл испытывает меньшую температурную деформацию, а также может содержать меньшую плотность дислокаций [107].
Обычно гидротермальная система состоит из высокого
цилиндрического сосуда, частично заполненного пористым питающим
веществом, растворителя, затравки (затравочный кристалл, выполняющий
роль центра кристаллизации) и перегородки. Мелкие частицы питающего
вещества находятся в нижней части сосуда, а затравки размещаются в его
верхней части. В сосуд помещается перфорированная металлическая
перегородка для разделения областей растворения и роста, и он
наполняется определенным количеством жидкости. Затем его вертикально
помещают в автоклав (аппарат для проведения различных процессов при
нагреве и под давлением выше атмосферного), разогреваемый таким
образом, чтобы нижняя часть с растворяющимся веществом была более
горячей, чем верхняя область роста кристалла. При определенном
давлении и температуре частицы вещества в нижней части сосуда
начинают растворяться, постепенно насыщая раствор, который
переносится конвективными течениями в зону роста. Кристалл растет в результате осаждения растворенного материала на пластинки затравок, так как вследствие более низкой температуры концентрация насыщения в этой области ниже. Процесс является медленным и для получения полноразмерного кристалла требуется несколько недель. На сегодняшний день существующие знания о росте кристаллов в гидротермальной системе являются результатом, полученным после многочисленных проб и ошибок [109].
В последнее время были предприняты попытки численного моделирования процессов, происходящих в автоклаве. Например, Roux и др. [117] были разработаны двумерная и трехмерная модели гидротермального роста. Однако в этой работе рассматривается сосуд с квадратным сечением вместо используемого в реальных системах цилиндрического автоклава. Chen и др. [108-110] была предложена комплексная 3-хмерная модель гидротермального роста. Конвективная система в автоклаве рассматривается как комбинация жидкого и пористого слоев. Однако расчеты проведены при отсутствии затравок и без учета изменения концентрации вещества в сосуде. В работе Li и др. [112] рассматриваются двумерные и трехмерные осесимметричные модели для изучения процессов в промышленных автоклавах, которые имеют большое соотношение высоты к диаметру - до 20. В результате трехмерного моделирования авторы заключили, что поток не является трехмерным, поэтому в статье представлены только результаты двумерного моделирования, В работе исследован перенос тепла через перегородку, сделан вывод, что он пренебрежимо мал, а перенос тепла практически полностью осуществляется потоком. Перенос массы, как и у предыдущих авторов, не учитывается.
В третьей главе настоящей работы рассматриваются процессы в
системе, которая используется для выращивания кристаллов берилла. До
сих пор нет полного представления о характере циркулирующих потоков в
сосуде, который подогревается со стороны боковой вертикальной стенки, и почти изолирован сверху и снизу. Ввиду большой сложности, и даже невозможности, провести измерения в реальной гидротермальной системе, температурное и концентрационное распределения остаются неизвестными. Единственная доступная информация - это размер и качество полученных кристаллов, а также затраченное на их рост время. Поэтому цель исследования на основе численной реализации модели [109,116] - рассмотреть конфигурации течения, температурного поля, а также изменение концентрации растворенного вещества выше пористой среды при гидротермальном росте кристаллов [100,103,115].
В первом параграфе настоящей главы описана математическая модель гидротермальной системы на основе объединения уравнений Навье-Стокса и Дарси-Бринкмана-Форчхеймера для получения общего решения в жидком и пористом слое. Так как процесс протекает очень медленно, то он предполагается квазистационарньш. Влияние растворимости и роста, а также концентрации питающего вещества на конвекцию не учитывается, Алгоритм численной реализации приведен во втором параграфе. В третьем параграфе представлены результаты анализа течения, распределения температуры и концентрации растворенного минерала.
Работа была выполнена при поддержке грантом РФФИ 02-05-64280.
Автор выражает искреннюю благодарность своему научному руководителю д.ф.-м.н. Попову Владимиру Николаевичу за постоянное внимание и руководство работой.
Результаты расчетов двумерных тестовых задач
Представим результат, полученный при реализации тестовой задачи о конвекции жидкости со свободной поверхностью в полости прямоугольного сечения. В недивергентной форме уравнения имеют вид
После предварительного преобразования исходных уравнений (1.13)-(1,15), была реализована следующая задача
Численное решение системы (1.19)-( 1.21) проводилось на равномерной пространственной сетке вида: x/=hx-i, i=0, ..., I, hx=\Il; zk=hz-k, h=0, ..., K, hz=]IK, которая разбивает расчетную область на IxK ячеек. По аналогии с методами типа MAC и SIMPLE [52,114], составляющие скоростей и, w определялись в серединах боковых граней ячеек. Давление Pi-]?2,k-\/2 рассчитывалось в центрах ячеек {x-hJ2, zk-hJ2\ /=1, .., І, к=\, К, а составляющие скоростей: u-Lk+\l2 в точках с координатами (х,-, zrhJ2), /=0, ..., I, =1, ..., К, _1/2д - в точках {X[-hJ2, zk), /=1, ..., /, =1,...Д(рис. 1.3).
Распределение давления и составляющих скоростей.
После аппроксимации балансных соотношений, получаемых интегрированием (1.19), (1.20), была построена система неявных разностных уравнений с порядком точности 0{x,hx,hz). Значения скоростей при определении 7/ н , 7/lif использовались с предшествующего временного слоя. Расчеты велись с двойной точностью, с временным шагом т=10"3 и прекращались после выхода рассматриваемого процесса на стационарный режим (Ы),5). На каждом шаге по времени для определения составляющих поля скоростей и давления применялся итерационный алгоритм типа Узавы [41,44,53]. Решение системы прекращалось на некоторой 5-ой итерации при одновременном выполнении условий max{«J-uJ"! } 10"s, тах{м/ -№ы } ]()Л i,k i,k
Результаты расчетов представлены в табл. 2, в столбцах которой приводятся величины максимальных отклонений точных от рассчитанных значений скоростей при Re=l и различных /, К. Отметим, что согласно полученных данных погрешность при использовании разностных операторов вида (1.5) для аппроксимации (1.19), (1.20) изменяется строго в зависимости h . Таблица 1хК тах{\ ulk-u(xhzk)\} і,к maM\wlk-w(xhzk)\} і,к 20x20 4.3380-10 3 8.8105-10 40x40 1.0529-10"3 2.0379-10"4 80x80 2.5928-10"4 4.8280-10"4
Представим результат, полученный с помощью описанного выше алгоритма для широко известной тестовой задачи о конвективном переносе тепла в полости прямоугольного сечения при подогреве сбоку (1.22) (1.23) (1.24) (1.25) [22,37] ди Зл 3iw ЗР 3 и д и 3 ск Зк + = + —- + & Зс Зс дР 3 w 3 w л „ + —- + —— + GrO, 3v 3uw 3w2 + З Зс 3: & де2 ck1 — + —= 0, Зс 3: -), 0 х 1, 0 z l, дв Зав 3v9 і у32в 32в З 3: ск Рг Зг А #(0,z,0 = 0, 0(W) = 1, M( W) = «(U,0 = 0 w(p,z,t) = w(l,z,t) = 0, 0 z l, 0(x,O,t) - x, &(x,l,t) = x, u(x,0,t) = u(x \ t) = 0, yv(x,0,t) = w(xXt) = 0, 0 х 1, ґ 0, w(x,z,0) = w(x,z,0) = 0, P(x,z,Q) = 0, в(х,г,0) = х, 0 x \, 0 z l. Здесь и, w - компоненты вектора скорости, Р - давление, 9 температура, Gr - число Грасгофа, Рг - число Прандтля. При определении решения предварительно было проведено преобразование исходных уравнений (] .22), (1.23), (1.25) к виду а ж ах ци ах, 3 т]ле az = - + 1(± ) + 1(± ) + Ог ?, (1.27) 3 dz Зс rjtl Зс az 7ДГ, az f [f( ) + f( )L (1.28) г -і - j -[rf -JWj -Pr \ud -Pr JW 7?„ - Є , 77и, = Є , //„ = Є , jUv. = Є, . Вычисление давления является частью процедуры решения задачи. Согласно [24], из (1.22)-(1.24) было получено уравнение ґдгР 32Р, Л2 nd2uw 32w2 Зв 3D -(—Г + т)=—Г + 2— + Г2—Gr—+ —, (1.29) ас 3: Зс 3x3 az az at _ ди Зл где D = -— + —. Зс dz Граничные условия для Р определяются с использованием (1.22) (1.24). Так, при 2=0: =_ + Ог0, (І.30) 3z 33с (аналогичное соотношение используется при z=l), а на боковых поверхностях при х=0 и х-1: = - -. (1.31) 3 33 Таким образом, система (1.29)-(1.31) представляет собой задачу Неймана для уравнения Пуассона, которую необходимо реализовьшать на каждом временном шаге. Решение (1.29)-(1.31) находилось с использованием преобразования Фурье для сеточных функций [35].
Использование экспоненциальных схем для решения трехмерных задач, описываемых уравнениями Навье-Стокса
Развитие многих отраслей современной науки и техники в значительной степени обусловлено получением и последующим использованием монокристаллов заданного состава, необходимого размера и высокого качества. Процессы конвективного тепломассопереноса в ростовых установках во многом определяют возможность получения кристаллов с заданными свойствами, поэтому поиску способов управления ими во всем мире уделяется большое внимание.
В условиях реальной технологии сложно или даже невозможно проводить достаточно подробные исследования взаимосвязанных процессов конвективного тепломассообмена в расплавах. Поэтому необходимо построение трехмерных математических моделей (так течения пространственные), алгоритмов их реализации, и проведение вычислительных экспериментов в широком диапазоне меняющихся параметров.
Несколько лет назад А.Е. Кохом была предложена, численно и экспериментально обоснована оригинальная идея выращивания кристаллов в азимутально распределенных стационарных и вращающихся (Heat Field Rotation Method - HFRM) тепловых полях. Суть этих методов состоит в том, что ростовой тигель разогревается вертикальными нагревательными элементами, равномерно расположенными по периметру его боковой стенки [72,73,90,91]. Нагревательные элементы включаются (в стационарном случае) или переключаются (HFRM) в некотором порядке, что изменяет температурное поле в жидкости и способствует возникновению в ней азимутальных течений, усиливая перемешивание. Это значительно упрощает технологию, так как отпадает необходимость переоснащения ростовой установки в разогретом состоянии - замена мешалки на з атрав ко держатель. Эффективность такого подхода демонстрируют результаты, полученные при выращивании кристаллов CsLiB6O!0 в циклически меняющемся тепловом поле [90].
В работе [80] рассмотрена структура течений в расплаве германия при циклически меняющемся неоднородном разогреве боковых стенок цилиндрического тигля. Результаты исследования подтвердили возможность применения бесконтактного метода для гомогенизации раствор-расплава. В [74,92,93] приводятся результаты численного моделирования нестационарной термогравитационной конвекции при получении монокристалла методом Чохральского в условиях циклически меняющегося температурного поля в расплавах кремния и германия. На основе анализа данных, полученных в ходе расчетов, сделаны выводы о существовании режимов теплового поля, обеспечивающих максимальное перемешивание кристаллизационной среды при Рг«1. Влияние HFRM на процессы тепломассопереноса в расплавах может способствовать разрушению диффузионного слоя у поверхности кристалла и уменьшению неоднородности распределения примеси в получаемых изделиях.
В настоящей главе рассматриваются математические модели процессов получения монокристалла методом Чохральского в азимутально распределенных стационарных тепловых полях и с применением HFRM, алгоритм их решения и некоторые результаты, полученные при исследовании зависимости конвективных течений в кристаллизационной среде обладающей свойствами Bi SiC o (BSO) с высокой вязкостью (Рг=26) от параметров теплового поля. Численное моделирование процессов конвекции во время роста кристалла основано на совместном решении нестационарных трехмерных уравнений Навье-Стокса в приближении Буссинеска и теплопереноса при заранее определенных условиях нагрева на границе расчетной области. Подробно описывается алгоритм реализации задачи с использованием экспоненциальной схемы.
2.1. Математическая модель ростовой установки для получения кристаллов методом Чохральского при неоднородном разогреве боковых стенок тигля
Схему рассматриваемой установки для выращивания кристаллов методом Чохральского иллюстрирует рис. 2.1. Расплавленный материал заполняет цилиндрический сосуд (тигель) с внутренним радиусом R0 до уровня Но. Из него вытягивается невращающиися монокристалл, радиус которого Rcr- Тигель помещен в нагревательную печь, средняя температура внутри которой То, с вмонтированными в ее боковые стенки поочередно включаемыми вертикальными нагревательными элементами. Разогрев жидкости происходит в результате теплообмена между стенками сосуда и окружающей средой, а отвод тепла осуществляется через ее свободную поверхность. Повышенный нагрев боковых стенок тигля в секторе площадью H0R0A(p определяется температурой Тн расположенного напротив включенного термоэлемента (Тн Т0).
Анализ температурного поля и структуры течений в условиях циклически изменяющегося разогрева боковых стенок ростового тигля
При проведении вычислительных экспериментов рассматривалась структура течений в раствор-расплаве Ві БіОго (ВSO) при дополнительном его разогреве от 12 поочередно включаемых вертикальных нагревательных элементов (НЭ), которые равномерно расположены вокруг сосуда (см. схему на рис. 2.1). Геометрические величины в расчетах использовались следующие: R j=Q.Q25 м, #0=0.03 м, RCr=0.0l м, ширина каждого нагреваемого сектора боковой стенки тигля лН0/6 (А р=7г/6). Физические свойства раствор-расплава BSO были определены согласно [95]: Х=0.345 Вт/(м-К), а=1.15-10"7 м2/с, р=7630 кг/м3, v=3-106 м2/с, р=7-10"5 1/К, 7 =1163 К. Значения температур - ГН 173 К, Гя=2003 К. В результате безразмерные параметры принимают значения: Pr=26, Gr=1.2-105, 0/-/=3, #=1.2, Rc=0A, Bi,=50, Ві2=2, Ві3=1.75.
Рис. 2.3-2.5 иллюстрируют результаты, полученные при рассмотрении варианта с несимметричным разогревом боковых стенок сосуда. На рис. 2.3д, 2.4а отражено квазистационарыое распределение температурного поля в жидкости при включенном НЭ с координатой фі=45. Видно, что максимальный разогрев раствор-расплава приходится на область около его свободной поверхности вблизи нагревательного элемента. Толщина этого слоя вдоль оси z незначительна, а перегрев не оказывает влияния на температурное поле в непосредственной близости от кристалла.
Проведем анализ структуры течений для рассматриваемого случая, Основная циркуляция в жидкости обусловлена вспльшанием теплых масс вверх вдоль боковых стенок, их охлаждением в приповерхностной области и формированием нисходящего потока вблизи центральной оси тигля. При этом наиболее интенсивные течения формируются непосредственно под кристаллом и у подогреваемого сектора боковой стенки (рис, 2.36). Ввиду того, что наибольший температурный градиент в жидкости в азимутальном направлении формируется в ее верхних слоях, именно здесь течение дополняется потоками, расходящимися от теплового источника в разные стороны. Значения азимутальной составляющей вектора скорости представлены на рис. 2.4с. Данные расчетов показывают существование азимутального течения в непосредственной близости от поверхности кристалла и быстрое снижение его интенсивности с увеличением глубины.
Первоначально, по результатам расчетов, было определено, что для влияния на структуру течений в жидкости со свойствами BSO ширина нагреваемого сектора на боковой стенке сосуда должна быть не уже %R0 /6, а температура теплового источника не ниже 9 =3, то есть требуется существенное тепловое воздействие. При меньших значениях этих параметров значительный температурный градиент в азимутальном направлении не возникал, а результат, подобный представленному на рис. 2.3-2.5, не удавалось получить даже качественно. Анализ влияния временного шага переключения нагревательных элементов Д на конвективные процессы в раствор-расплаве показал, что устойчивый поток в азимутальном направлении формируется при значениях этого параметра не меньше 2. Рис. 2.5а,с иллюстрируют траектории движения маркеров в жидкости, а рис. 2.5b,d- их проекции на плоскость гц при переключении нагревательных элементов, включаемых через временной интервал Д=2 в направлении против часовой стрелки. Согласно представленным результатам, происходит устойчивое азимутальное перемещение части жидкости в направлении, совпадающем с направлением, в котором переключаются тепловые источники. Этот эффект обусловлен тем, что существует синхронизация режима переключения нагревательных элементов с периодом циркуляции раствор-расплава между боковыми стенками и центральной частью сосуда. а Ь
Наблюдаемые частицы в то время, когда они находятся в придонной области, переносятся потоком строго по прямой от центра сосуда к его стенке. Далее, поднимаясь вверх, они начинают перемещаться к центру тигля. Приобретаемая при этом в зоне перегрева величина азимутальной составляющей вектора скорости сильно зависит от расстояния маркера до поверхности жидкости. Так, если рассматриваемая частица попадает в область, далекую от стенок тигля, и не достигает верхних слоев раствор-расплава, ее движение осуществляется почти строго в плоскости rz (более плотные участки траектории на рис. 2.5 з, 2.56). В случаях, когда рассматриваемый объем жидкости поднимается вверх в непосредственной близости от боковых стенок, он движется к центру сосуда почти у самой поверхности. Азимутальное отклонение здесь максимальное
Численные исследования процессов гидротермального синтеза при получении искусственных изумрудов
Существует множество параметров, влияющих на течение в автоклаве, например, его размер, распределение температуры вдоль стенок, высота пористого слоя, глубина, на которую опускается пластина затравки, и ее ширина, число Дарси, число Грасгофа и т. д. В настоящем параграфе представлены результаты численных решений, которые позволяют проиллюстрировать возможные структуры потока и полей температуры и концентрации в сосуде.
При проведении вычислительных экспериментов рассматривались структура течений и температурные поля в растворе и пористой среде. В виду отсутствия четкой информации о теплофизических свойствах растворяющегося вещества и растворов в условиях высокого давления при получении искусственных изумрудов, при расчетах использовались следующие параметры. Физические свойства водного раствора при высоком давлении были определены согласно [98,99,101]: ц„=0.54 Вт/(м-К), с =5736 Дж/(кг К), a =1.3240-7 м2/с, pw=712.5 кг/м3, v=1.28-10 7 м2/с, [3=2.92-10-3 1/К. Для растворяющегося вещества использовались теплофизические свойства берилла (Ве3А128ібОі8) р/=2650 Дж/(кг-К), 0,=1255 Дж/(кг К), Хь=6.276 Вт/(м К) (аь=Хь/(ръ сь)) [11.8]. Скорость роста кристалла w0 была оценена как 10"4 м/сутки [96,119]. Распределение температуры вдоль боковой стенки сосуда определено в результате замеров на работающей установке для получения изумрудов и представлено на рис. 3.2. Максимальная Ть и минимальная Tt температуры вдоль вертикальной стенки сосуда 903 К и 853 К соответственно.
Геометрические величины в расчетах согласно измерениям на рабочей установке: /?о=0.0.15 м, высота сосуда - 0.3 м, уровень, выше которого располагается затравка - 0.15 м, высота пористого слоя - 0.12-0.15 м, ширина пластины затравки - 0.01875 м. Величины безразмерных параметров, используемых при расчетах, были определены как: #=20, Sc=128, Pr=l, Gr=106, а=10-3-10"6, Ві0=1.36-Ш4, &=є+1.23(1-є), y?t=-rl 4.94(1-є). Значение пористости "-0.4, а соотношение вязкостен Л = 1/е = 2.5 [109,110]. Высота пористого слоя #;=8-10, координата, выше которой располагалась пластина затравки, Д=10, а отношение ширины затравки к внутреннему диаметру сосуда составляло 0,625.
Вычисления проводились на пространственных сетках со значениями I JxK от 24x24x90 до 40x72x120 с хорошим качественным и количественным совпадением получаемых результатов. Значения
временного шага выбирались из условий устойчивости численного счета, минимизации количества итераций и изменялись от 10 J до 2-Ю"5.
Так как гидротермальный процесс выращивания кристаллов является медленным, то его на некотором временном интервале можно считать квазистационарным и анализировать поля температуры и потока для заданной высоты пористого слоя, пренебрегая влиянием растворимости и роста. Рассмотрим конфигурацию течения, температурного поля, а также изменение концентрации растворенного вещества выше пористой среды для некоторого момента времени. Рис. 3.3, 3.4 отображают поле скоростей (а), изотермы (Ь) и изоконцентраты (с) при Нр=\0 в плоскостях rz при ф=2.5, 45 к поверхности затравки, а также в поперечном сечении автоклава.
Течение в жидком слое является многоячеистым и пространственным. В области 10 z 17 имеется несколько завихрений. В верхней части при z 17, где температурный градиент вдоль вертикальной оси z почти отсутствует, скорость течения существенно снижается. Стрелки на рис. 3.3а и 3.4Й указывают направление, а их длина характеризует интенсивность потока. Прямая линия в основании полукруга на рис, 3.4 проходит через плоскость затравки.
Циркуляция жидкости выше пористого слоя (рис. 3.3а) обуславливается наличием завихрения вблизи пористого слоя, которое формируется всплыванием быстро нагревающегося раствора вдоль боковых стенок (/=1). При z«12 поднимающийся поток сталкивается со спускающимся сверху и заставляет его отклониться к центральной части сосуда, где формируются нисходящее и восходящее течения. Наибольшую интенсивность потоки имеют на некотором удалении от поверхностей боковых стенок сосуда и пластины затравки.
