М ЕХАНИКА ЖИДКОСТИ И ГАЗА № 3 • 2015
УДК 532.511:536.252
КОНВЕКТИВНЫЕ ВЗАИМОДЕЙСТВИЯ И УСТОЙЧИВОСТЬ ТЕЧЕНИЙ В МОДЕЛИ МЕТОДА ЧОХРАЛЬСКОГО ПРИ ВРАЩЕНИИ КРИСТАЛЛА
© 2015 г. О. А. БЕССОНОВ
Институт проблем механики им. А.Ю. Ишлинского РАН, Москва e-mail: bess@ipmnet.ru
Поступила в редакцию 22.10.2014 г.
Представлены результаты расчетов по исследованию влияния вращения кристалла на устойчивость течений в широком диапазоне значений числа Прандля (Pr). Для низких Рг определены режимы с повышенным порогом устойчивости, рассмотрены механизмы потери устойчивости при превышении критических значений числа Грасгофа (Gr) и скорости вращения кристалла. Для средних и высоких Рг представлены таблицы режимов течений, рассмотрены особенности устойчивых неосесимметричных спиральных течений различных типов, выявлены зоны частичной стабилизации течения.
Ключевые слова: тепловая гравитационная конвекция, рост кристаллов из расплава, гидродинамическая модель метода Чохральского, численное моделирование, конвективная устойчивость, конвективные взаимодействия, стабилизация колебаний.
Метод Чохральского для выращивания монокристаллов из расплавов характеризуется широким набором параметров, влияющих на устойчивость течения и, как следствие, на качество получаемых кристаллов. К их числу относятся физические характеристики жидкости, геометрические параметры установки, тепловые граничные условия, параметры свободной и вынужденной конвекции. Исследование условий и механизмов потери устойчивости течения при комплексном влиянии различных факторов позволяет определить оптимальные режимы и устранить либо ослабить колебания температуры в расплаве, что необходимо для снижения уровня неоднородностей распределения примесей в выращиваемых кристаллах.
Для изучения различных механизмов течения и их взаимодействия широко используется гидродинамическая модель метода Чохральского, характеризующаяся упрощенной геометрией и граничными условиями, что позволяет использовать эффективные численные методы и проводить расчеты с хорошим пространственным и временным разрешением [1]. Для сравнения численных методов и вычислительных кодов ставятся тестовые задачи на основе этой модели, что повышает надежность численного исследования [2, 3]. Упрощенная модель установки Чохральского также широко используется в экспериментальных исследованиях [4—6].
Разработанные в последние годы вычислительные коды, а также развитие вычислительной техники позволили проводить моделирование сложных пространственных нестационарных течений. Однако детальное изучение конвективных взаимодействий затрудняется наличием большого числа параметров задачи. По этой причине систематические исследования проводились, как правило, в постановке с неподвижными тиглем и кристаллом [1, 7—9]. Такое упрощение позволило провести масштабное параметрическое исследование конвективных взаимодействий в широком диапазоне чисел Прандтля, выявить основные механизмы течения и построить карту режимов для
Сгс
108
107
106
0.
Фиг. 1. Зависимость критического числа Грасгофа Огс от Рг для осесимметричной (1) и пространственной (2) моделей; характеристики течения: А — стационарное в обеих моделях, В — нестационарное в пространственной модели, С — всегда нестационарное, Б — нестационарное с сохранением осевой симметрии в пространственной модели
осесимметричного и пространственного режимов моделирования [10, 11]. Исследования с вращением кристалла и(или) тигля также проводились, но, как правило, изучались конкретные режимы в узких диапазонах параметров [5, 6, 11—13].
Настоящая статья продолжает цикл исследований, выполненных в предыдущие годы под руководством проф. В.И. Полежаева (1936—2013 гг.). В статье приводятся результаты параметрического исследования, проведенного на усложненной модели с учетом вращения кристалла при неподвижном тигле в широком диапазоне чисел Прандтля. Выбор параметров моделирования производится на основе полученной ранее обобщенной карты режимов тепловой гравитационной конвекции [10, 11]. Основная цель работы — нахождение режимов с повышенным порогом устойчивости и зон стабилизации течения, а также рассмотрение механизмов взаимодействия тепловой гравитационной и вынужденной конвекций с построением сводных таблиц режимов.
1. Математическая модель и численный метод. Установка гидродинамической модели метода Чохральского представляет из себя цилиндрический тигель радиуса Яс, заполненный расплавом до высоты Н (фиг. 1, [9]). В центре к поверхности расплава примыкает диск радиуса имитирующий торец выращиваемого кристалла. Открытая поверхность расплава предполагается плоской.
Численное моделирование основано на решении пространственных нестационарных уравнений Навье—Стокса в приближении Буссинеска для несжимаемой вязкой жидкости в переменных скорость V, давление р, температура 9
ддУ + V • (уу) = -Vр + V 2У - в^е д г
V-у = 0 (1.1)
+ (У0) = — V 20 д г Рг
Здесь вг = gв ЬМ/V2 — число Грасгофа, Рг = v/a — число Прандтля, L и AT — масштабы длины и температуры, в — коэффициент температурного расширения, V — кинематическая вязкость, а — температуропроводность, g — вектор силы тяжести. Масштаб скорости определяется как У0 = у/Ь, масштаб времени — г0 = Ь /V. Используется цилиндрическая система координат 0 <ф< 2п, 0 < z < Н, 0 < г < Яс.
В соответствии с постановкой задачи в [10] геометрические параметры расчета определены следующим образом: Rc = 1, H = 1, Rx = 0.4. На поверхности расплава выполняются условия скольжения: дУг/дz = 0, д У^/дz = 0, Уг = 0. Кристалл вращается с угловой скоростью Яе^. Используются следующие граничные условия для температуры: нагрев на боковой цилиндрической границе 0 = 1, охлаждение на границе расплава и кристалла 0 = 0, тепловая изоляция на дне тигля и на поверхности расплава 59/д1 = 0.
Для решения уравнений (1.1) используется метод конечных объемов. Применяются разнесенные сетки со сгущениями в осевом направлении. Дискретизация конвективных членов производится по схеме с квадратичными интерполяциями против потока. Для температуры 9 при высоких Рг, а также для переноса азимутальной компоненты скорости У^ в радиальном направлении применена монотонная неосциллирующая схема с ограничителем потока. Уравнения решаются раздельно, с использованием метода проекций и частично неявной схемы интегрирования по времени. Для решения уравнения Пуассона для давления применен экономичный прямой метод с разложением по базису собственных функций оператора Лапласа [14]. Расчетная программа распараллелена для многопроцессорных и многоядерных вычислительных систем с общей памятью в модели ОрепМР.
Численный метод является прямым и не содержит дорогостоящих итерационных этапов. Благодаря высокой эффективности он может использоваться для детального моделирования ламинарных, переходных и турбулентных пространственных течений на подробных сетках. Результаты сравнения с экспериментом показывают, что используемый метод обеспечивает высокую точность и реалистичность моделирования сложных динамических пространственных течений. Развернутое изложение метода представлено в [14].
Все расчеты производились в пространственном режиме моделирования. Использовались сетки размером 64 х 87 х 80 и 64 х 120 х 120 (ф, z, г). Для высоких значений вг расчеты уточнялись на сетке 128 х 120 х 120. Сетки имеют сгущение в пропорции 1:4 около дна полости и у поверхности расплава.
2. Постановка задачи моделирования. Тепловая гравитационная конвекция расплава в тигле вызывается действием двух механизмов: нарушением механического равновесия при подводе тепла к боковой поверхности тигля и потерей устойчивости равновесия в связи с отводом тепла от торца кристалла. Первый механизм, характерный для низких Рг, вызывает глобальное или подъемно-опускное течение, захватывающее весь объем жидкости. Второй механизм начинает действовать при повышении Рг, что приводит к снижению теплового потока на границе кристалла и возникновению неустойчивости типа Рэлея—Бенара. Этот механизм вызывает локальное течение, имеющее
характер струй либо ячеек (термиков), опускающихся от границы кристалла вниз. Вынужденная конвекция, возникающая при вращении кристалла, вступает во взаимодействие с этими механизмами.
Главная цель настоящего исследования — определение сочетаний значений числа Грасгофа и скорости вращения кристалла (Ог, Яе^), при которых происходит потеря устойчивости либо смена режима течения, а также изучение режимов течений в широком диапазоне чисел Прандтля. Вращение кристалла рассматривается в первую очередь как способ полной либо частичной стабилизации течения. Необходимость вращения кристалла может быть также вызвана технологическими причинами, поэтому представляет интерес нахождение предельных значений Яех, при которых сохраняется стационарный режим.
В качестве характеристики взаимодействия вращения кристалла и тепловой конвекции используется параметр Яе^/Ог, определяющий соотношение интенсивностей этих механизмов. При значении параметра порядка единицы механизмы вынужденной и тепловой конвекций частично компенсируют друг друга. С учетом соотношения радиусов кристалла и тигля пороговое значение параметра взаимодействия для рассматриваемой конфигурации находится в интервале от 2 до 2.5 [11]. Как правило, при превышении этого порога устойчивость течения снижается, поэтому ниже рассматриваются умеренные значения параметра Яе^/Ог, при которых вклад вращения кристалла ниже вклада тепловой конвекции или сопоставим с ним.
На фиг. 1 приведена сводная карта режимов неустойчивости тепловой гравитационной конвекции в виде диаграммы (Рг, Огс) для варианта модели метода Чохральского с теплоизолированной поверхностью расплава в отсутствие вращений кристалла и тигля [10, 11].
Для изучения влияния вращения кристалла на устойчивость течения на диаграмме выбраны 4 позиции, характеризующие различные типы и режимы течений: Рг = 0.03 (а) — типичное течение с глобальным механизмом конвекции; Рг = 0.15 (б) — точка стабилизации, соответствующая максимальному критическому числу Огс; Рг = 1 (в) — смена моды потери устойчивости течения с трехмерной
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.