CC BY
9
i L Научно-технические ведомости СПбГПУ. Физико-математические науки. 15 (3) 2022
St. Petersburg State Polytechnical University Journal. Physics and Mathematics. 2022 Vol. 15, No. 3 --►
Математическое моделирование физических процессов
Научная статья УДК 532.5
DOI: https://doi.org/10.18721/JPM.15303
СРАВНИТЕЛЬНАЯ ОЦЕНКА RANS-МОДЕЛЕЙ
ТУРБУЛЕНТНОСТИ С ИЗОТРОПНОЙ ВЯЗКОСТЬЮ
ДЛЯ РАСЧЕТА КОНВЕКЦИИ РАСПЛАВА КРЕМНИЯ
В УСТАНОВКАХ ВЫРАЩИВАНИЯ КРИСТАЛЛОВ
Д. В. Борисов н, В. В. Калаев
АО «Группа СТР» - ООО «Софт-Импакт», Санкт-Петербург, Россия
н dmitriy.borisov@str-soft.com
Аннотация. В работе сопоставляются результаты RAN S-расчетов турбулентной конвекции в расплаве кремния, полученные по нескольким моделям турбулентности с изотропной вязкостью, с ранее опубликованными данными вихреразрешающих ILES-вычислений для аналогичных условий. Выбирается модель турбулентности для ее последующей проблемно-ориентированной модификации с алгебраическим введением факторов, которые могут продуцировать нужную анизотропию тензора рейнольдсовых напряжений и вектора турбулентного теплового потока, входящих в осредненные по Рейнольдсу уравнения движения и энергии. Показано, что применительно к задачам расчета конвекции в тиглях установок, где используют метод Чохральского, целесообразно взять для модификации либо однопараметрическую k-модель, либо двухпараметрическую k-s модель в качестве исходной RANS-модели.
Ключевые слова: метод Чохральского, ILES, RANS, тензор рейнольдсовых напряжений, турбулентный теплоперенос, конвекция
Для цитирования: Борисов Д. В., Калаев В. В. Сравнительная оценка RANS-моделей турбулентности с изотропной вязкостью для расчета конвекции расплава кремния в установках выращивания кристаллов // Научно-технические ведомости СПбГПУ. Физико-математические науки. 2022. Т. 15. № 3. С. 28-42. DOI: https://doi.org/10.18721/ JPM.15303
Статья открытого доступа, распространяемая по лицензии CC BY-NC 4.0 (https:// creativecommons.org/licenses/by-nc/4.0/)
Original article
DOI: https://doi.org/10.18721/JPM.15303
COMPARATIVE EVALUATION OF RANS EDDY-VISCOSITY TURBULENCE MODELS FOR CALCULATING THE SILICON MELT CONVECTION IN CRYSTAL GROWTH SYSTEMS D. V. Borisov , V. V. Kalaev
STR Group, Inc. - Soft-Impact, Ltd., St. Petersburg, Russia
H dmitriy.borisov@str-soft.com
Abstract. In the paper, the results of RANS calculations of turbulent convection in silicon melt, obtained using several eddy-viscosity turbulence models, have been compared with previously
© Борисов Д. В., Калаев В. В., 2022. Издатель: Санкт-Петербургский политехнический университет Петра Великого.
published ILES eddy-resolving calculation data for similar conditions. A turbulence model was chosen for its subsequent problem-oriented modification with the algebraic introduction of factors that could produce the required anisotropy of the Reynolds stress tensor and the turbulent heat flux vector including in the Reynolds-averaged equations of momentum and energy. As applied to the problems of calculating the convection in crystal growth furnace crucibles using the Czochralski method, it was shown the expediency of taking either the one-equation k-model or the two-equation k-s model as the initial RANS-model for modification.
Keywords: Czochralski method, ILES, RANS, Reynolds stress tensor, turbulent heat transfer
For citation: Borisov D. V., Kalaev V. V., Comparative evaluation of RANS eddy-viscosity turbulence models for calculating the silicon melt convection in crystal growth systems, St. Petersburg State Polytechnical University Journal. Physics and Mathematics. 15 (3) (2022) 28-42. DOI: https://doi.org/10.18721/JPM.15303
This is an open access article under the CC BY-NC 4.0 license (https://creativecommons. org/licenses/by-nc/4.0/)
Введение
Метод Чохральского является одним из основных методов производства полупроводниковых кристаллов кремния, широко применяемых в электронной промышленности [1 — 3]. Для получения монокристаллов высокого качества необходимо контролировать массоперенос примесей в расплаве, попадающих в кристалл во время ростового процесса. Течение расплава кремния, как правило, — турбулентное, что задается условиями выращивания кристаллов даже относительно малого диаметра (10 см) в лабораторных установках. Наличие разномасштабных турбулентных структур в тиглях для промышленного выращивания кристаллов диаметром 20 —30 см затрудняет контроль концентрации примесей, а также может приводить к смене монокристаллического режима роста поликристаллическим [4]. Экспериментальное исследование турбулентного течения расплава кремния затруднено как высокими температурами процессов, так и требованиями к точности оборудования, используемого для измерения турбулентных пульсаций. В связи с этим, на сегодняшний день численное моделирование представляется наиболее перспективным методом исследования турбулентного течения и процессов тепло- и массообмена в расплаве кремния.
Наиболее точным методом расчета турбулентных течений является прямое численное моделирование (англ. Direct Numerical Simulation (DNS)), призванное разрешать все пространственно-временные масштабы турбулентности и не требующее дополнительных гипотез для замыкания уравнений [2, 3]. Однако данный метод требует значительных вычислительных ресурсов, что делает невозможным его применение в практических инженерных расчетах. Наиболее востребованным и относительно экономичным на сегодняшний день является подход, использующий уравнения Навье — Стокса, осредненные по Рейнольдсу (англ. Reynolds-averaged Navier — Stokes (RANS)), и позволяющий проводить расчеты в осесимметричной постановке (в том числе).
Для того чтобы численно изучать влияние вращения тигля на теплообмен в расплаве кремния на основе RANS-подхода, авторы работы [5] использовали k-s модель турбулентности. Результаты расчетов показали общее качественное согласие с экспериментальными данными. В работе [6] выполнены расчеты с использованием низкорейнольдсо-вой k-s модели, которые выявили, в частности, факт усиления тепло- и массопереноса по мере увеличения скорости вращения тигля. Подобных результатов не удавалось достичь при проведении предыдущих расчетов, заданных в ламинарной постановке, которая не подразумевает использования какой-либо модели турбулентности. Было показано также хорошее согласие расчетных распределений температуры в расплаве и формы интерфейса с экспериментальными данными. Авторы статьи [7] провели тестирование трех моделей турбулентности для расчета турбулентной конвекции расплава: «стандартной» k-s модели с использованием пристенных функций; двухслойной k-s модели, сочетаемой с
© Borisov D. V. Kalaev V. V., 2022. Published by Peter the Great St.Petersburg Polytechnic University.
моделью одного уравнения у твердых стенок; и низкорейнольдсовой k-s модели Лаундера — Джонса [8]. Было продемонстрировано преимущество последней из указанных моделей, которое заключалось в способности давать решения, близкие к ламинарным, в случае вырождения турбулентности. В работе [9] были исследованы турбулентные характеристики расплава кремния в идеализированном цилиндрическом тигле с применением нестационарного RANS-подхода (Unsteady RANS (URANS)) в трехмерной постановке. В своих расчетах они использовали k-s модель Лаундера — Шармы [10], а также k-ю SST-модель Ментера [11, 12]. Результаты этих URANS-расчетов показали преимущество SST-модели, которое заключается в лучшем разрешении структуры течения, температурных градиентов в окрестности стенок, а также получении более интенсивного течения на свободной поверхности расплава.
При численном моделировании конвекции расплава кремния в тиглях установок (работа по методу Чохральского), проведенном на основе осесимметричной RANS-формули-ровки, определенные успехи были достигнуты. Однако ряда важных характеристик тепло-и массопереноса в расплаве воспроизвести в рамках этого подхода не удалось. Например, для корректного моделирования термических напряжений и точечных дефектов в объеме кристалла необходимо правильно предсказывать форму фронта кристаллизации, которая существенно зависит от структуры течения и теплообмена в расплаве. В частности, в работах [13 — 16], где для моделирования теплопереноса в расплаве использовали RANS-подход, говорится о различии в 2 — 3 раза экспериментального и расчетного прогибов интерфейса. Ключевая причина такого рассогласования для кристаллов диаметром 100 мм заключается в том, что RANS-расчеты предсказывают интенсивное нисходящее течение в окрестности оси симметрии тигля, которое не наблюдается экспериментально [17]. Авторы работы [16] вводили специальные поправки в RANS-модель турбулентности для улучшенного предсказания формы фронта кристаллизации при моделировании процесса роста кристаллов диаметром 300 мм, создаваемых методом Чохральского. Однако, как отмечают сами авторы, применение данных поправок было направлено на частное, искусственное устранение недостатка RANS-модели турбулентности в подкристальной зоне, а не на общее улучшение собственно модели.
Другой важной характеристикой конвекции в расплаве кремния, влияющей на свойства и качество кристалла, выступает концентрация кислорода в расплаве. Предпринятые ранее попытки использования RANS-модели турбулентности для предсказания уровня концентрации кислорода в широком диапазоне параметров, управляющих работой установки по методу Чохральского, не приводили к результатам, согласующимся с опытными данными [6, 18].
Большинство RANS-моделей исходит из гипотезы Буссинеска, основанной на концепции изотропной турбулентной вязкости и из стандартной гипотезы градиентной диффузии (англ. Standard Gradient Diffusion Hypothesis (SGDH)), замыкающей осредненные по Рейнольдсу уравнения переноса температуры и/или концентрации.
В работе [19] мы исследовали вопросы локальной применимости гипотез Буссинеска и SGDH для моделирования тензора рейнольдсовых напряжений и вектора турбулентного теплового потока в расплаве кремния. Для этого были специально обработаны и затем проанализированы результаты вычислений, выполненных на основе вихреразрешающего метода (англ. Implicit LES (ILES)). При этом было установлено, что сильная анизотропия тензора рейнольдсовых напряжений в окрестности стенки тигля, интерфейса кристаллизации, а также у свободной поверхности расплава не воспроизводится на основе гипотезы Буссинеска; последняя исходно нацелена лишь на описание сдвиговых напряжений. Кроме того, выраженная анизотропия турбулентного теплопереноса у свободной поверхности также не описывается в рамках гипотезы SGDH.
Основные цели данного исследования формулируются следующим образом:
сопоставить результаты RANS-расчетов турбулентной конвекции в расплаве кремния, полученные по нескольким моделям турбулентности с изотропной вязкостью, с данными, полученными в рамках вихреразрешающего подхода ILES и опубликованными ранее в работе [19];
выбрать RANS-модели турбулентности для последующих их модификаций с настройкой на основе данных ILES-расчета.
Представляемые в данной статье результаты RANS-расчетов сопряженного теплообмена получены с применением программного пакета Flow Module, который является частью программного продукта CGSim, разработанного в STR Group и предназначенного для моделирования тепло- и массопереноса при выращивании кристаллов различными методами [20].
Мы использовали следующие низкорейнольдсовы модели турбулентности: k-модель Вольфштайна [21], k-s модель Чена [22] и k-ю SST-модель Ментера [11, 12].
Применяя пакеты Flow Module и ANSYS Fluent для случая модели Ментера, мы также провели кросс-верификационные расчеты в условиях упрощенной задачи: с моделированием турбулентной конвекции только расплава.
Математическая модель
При численном моделировании процесса выращивания кристаллов кремния методом Чохральского рассматривается турбулентная конвекция и теплообмен в расплаве и в инертном газе (аргон), который прокачивается над расплавом. Одновременно рассчитывается теплопроводность в кварцевом и графитовом тиглях, а также в монокристалле кремния. Математическая модель на основе RANS-подхода включает стационарные уравнения баланса массы, импульса и температуры:
V-(pu ) = 0, (1)
V- (рйй) = - V + V- (т - pu'u') + (р-Ро) g (2)
т = ц ( VU + ( VU) ) -1 ц ( V - U) E, (3)
V- (p cp ) = V- ( XVT - p cpT ), (4)
p (T ), если расплав или твердый блок;
Р = РоМ (5) ——, еслигаз.
RT
Здесь р, р кг/м3, — локальная и равновесная плотности, соответственно; u, м/с, — вектор скорости; р, р Па, — локальное давление и давление газа в ростовой установке, соответственно; g — ускорение свободного падения; т — тензор вязких напряжений; Е — единичный тензор; ц, Па-с, — динамическая вязкость; с , Дж/(кг-К), — удельная теплоемкость при постоянном давлении; Т, К, — температура; X, Вт/(м-К), — коэффициент теплопроводности; М, а.е.м., — молекулярная масса; R , Дж/(кмоль-К), — универсальная газовая постоянная.
Верхняя черта означает осреднение величины по Рейнольдсу, верхний штрих — пуль-сационную величину; и'и' и и'Т' означают тензор рейнольдсовых напряжений и вектор турбулентного теплового потока, соответственно.
Турбулентный перенос импульса и тепла моделируются в рамках концепции изотропной турбулентной вязкости и стандартной гипотезы градиентной диффузии:
2, ( - t ~\T 2t -\ \
u'u'=- kE- v 3
Vu + (Vu) —(V-u) E
(6)
и'Т' = - £ ^ (7)
1 ч
где м2/с, — кинематическая турбулентная вязкость, которая определяется в соответствии с моделью турбулентности; Рг^ — турбулентное число Прандтля.
Система уравнений (1) — (5) дополняется одним или двумя уравнениями переноса: кинетической энергии турбулентности k, м2/с2 (для всех использованных моделей), скорости диссипации 8, м2/с3 (для модели Чена) и удельной скорости диссипации ю, 1/с (для модели Ментера).
Согласно моделям Вольфштайна и Чена, турбулентная вязкость определяется следующим образом:
п гк2
V = С/,— > (8)
8
где С — эмпирическая константа, / — демпфирующая функция.
Согласно 88Т-модели Ментера, турбулентная вязкость определяется как
V, =-^-, (9)
тах(а1ш, SF2)
где a1 — константа модели; S, 1/с, — модуль тензора скоростей деформаций, F2 — функция модели, введенная в работе [11].
Вычислительные средства
Расчеты турбулентного теплопереноса в расплаве проводились в двумерной/осесим-метричной версии программного пакета Flow Module. Для дискретизации системы уравнений (1) — (5) применялся метод конечных объемов. Решение дискретизированных уравнений осуществлялось на основе алгоритма SIMPLEC. Расчет конвективных потоков температуры и компонент скорости выполнялся с использованием схемы QUICK, а кинетической энергии турбулентности и удельной диссипации — с использованием противо-поточной схемы первого порядка точности. Расчет диффузионных потоков осуществлялся с точностью второго порядка.
При проведении кросс-верификационных расчетов для представленной ниже упрощенной модельной задачи использовалась также двумерная/осесимметричная версия программного пакета ANSYS Fluent. Использовался алгоритм SIMPLEC для решения дискретизированных уравнений, расчет конвективных потоков компонент скорости и температуры проводился по схеме QUICK, конвективные потоки к и ю рассчитывались по противопоточной схеме первого порядка точности.
Модельная задача
Расчетная область модельной задачи включала в себя только область расплава (рис. 1). Радиус тигля R составлял 170 мм, радиус кристалла R = 50 мм, высота расплава H = 97 мм.
Рис. 1. Расчетные область (справа) и сетка (слева) модельной задачи:
АВ — граница раздела расплав — кристалл; ВС, CD — поверхности раздела расплав — газ (свободная) и расплав — тигель; ОА — ось симметрии; Rc, Rs — радиусы расплава и кристалла соответственно; Н — высота расплава; юс, ю^ — угловые скорости вращения тигля и кристалла соответственно
Ставились следующие граничные условия: на стенке тигля и фронте кристаллизации использовалось условие прилипания; на свободной поверхности расплава задавалось нулевое сдвиговое напряжение; на стенке тигля задавалась зависимость температуры от отношения радиуса к глубине расплава, полученная из экспериментальных значений [17]; фронт кристаллизации брался при постоянной температуре, равной температуре плавления кремния; свободная поверхность считалась адиабатической.
Свойства жидкого кремния, использованные в расчетах, представлены в табл. 1. Модельная задача решалась без учета эффекта Марангони. Скорости вращения стенки тигля юс и стенки кристалла ю^ составляли —5 и 20 об/мин соответственно (см. рис. 1).
Расчетная сетка содержала около 72 тыс. ячеек, размер первой приповерхностной ячейки составлял около 0,12 мм, размер ячейки в объеме — около 1,8 мм. Безразмерная координата у+ в первой пристенной ячейке не превышала единицы.
Таблица 1
Значения параметров жидкого кремния, использованные в расчетах
Параметр Обозначение Единица измерения Значение параметра
Плотность Р кг/м3 3194 - 0,3701 Т
Равновесная плотность Ро кг/м3 2570
Теплопроводность X Вт/(мК) 66,5
Теплоемкость с Дж/(кгК) 915
Динамическая вязкость ц Пас 8-10"4
Температура плавления T K 1685
Коэффициент Марангони до/дТ Н/(мК) -110-4
Степень черноты srad - 0,3
Примечание. Данные, приведенные в двух последних строках, будут использованы далее.
Рис. 2. Сравнения радиальных распределений температуры (а); осевой (b), радиальной (с), окружной (d) компонент скорости, а также турбулентной вязкости (e), рассчитанных при разной глубине расплава: 1 см (1, 3) и 3 см (2, 4) и при использовании программных пакетов
ANSYS Fluent (1, 2) и Flow Module (3, 4)
На рис. 2 представлены сравнения радиальных распределений температуры, компонент скорости и турбулентной вязкости на глубине расплава 1 и 3 см, полученных с использованием программных пакетов Flow Module и ANSYS Fluent. Видно, что результаты, полученные с использованием обоих пакетов, хорошо согласуются между собой. Можно отметить небольшие различия в распределениях осевой и радиальной компонент скорости в окрестности оси симметрии, что может быть обусловлено различными аппроксимациями слагаемых, обратно пропорциональных радиальной координате; слагаемые возникают при записи уравнения баланса импульса в цилиндрической системе координат.
100 г-
Angle, degrees
Рис. 3. Распределения плотности теплового потока вдоль тигля, полученные с помощью
ANSYS Fluent (1) и Flow Module (2)
На рис. 3 приведено сравнение распределений плотности теплового потока вдоль стенки тигля, полученных с помощью пакетов Flow Module и ANSYS Fluent. Разница между распределениями нарастает по мере приближения к оси симметрии, что может быть следствием различий в структуре течения. Несмотря на небольшое расхождение между результатами, полученными с использованием обоих кодов, можно сделать вывод о корректности имплементации модели Ментера в программный пакет Flow Module.
Постановка сопряженной задачи
Расчетная область сопряженной задачи, сформулированной на основе данных для ростовой установки EKZ-1300 [23], включает расплав, кристалл, кварцевый и графитовый тигли, а также часть газовой области над расплавом. Схема расчетной области приведена на рис. 4. Радиус тигля Rc составлял 170 мм, радиус кристалла Rs = 50 мм, высота расплава H = 97 мм.
Рис. 4. Схема расчетных области (справа) и сетки (слева) сопряженной задачи:
1 — расплав; 2 — кристалл; 3, 4 — кварцевый и графитовый тигли, соответственно; 5 — течение аргона; 6 — ось симметрии; 7 — свободная поверхность (остальные обозначения те же, что на рис. 1)
Ставятся следующие граничные условия: на границах раздела расплав-тигель и расплав-кристалл задается условие прилипания, на выходной границе — нулевое выходное давление.
Граничное условие на свободной поверхности расплава учитывает термокапиллярный эффект Марангони:
( ди Л ( ди Л дс дТ
к
дп
melt
к
дп
+ -
gas
дТ дт
(10)
т соответствует направлению, касательному к свободной поверхности, а п — нормальному; индексы melt и gas обозначают расплав и газ, соответственно.
Фронт кристаллизации поддерживается при постоянной температуре, равной температуре плавления кремния.
На внешних границах используется следующее граничное условие:
1 =
дп
дп
Т4 -Q
i ext
rad'
(11)
где gsb — константа Стефана — Больцмана; Emd — степень черноты; Q'"ad — падающий радиационный тепловой поток, полученный из решения задачи глобального теплообмена для установки EKZ-1300 (рис. 5); индекс ext соответствует внешней границе, индекс gas — соседнему газовому блоку.
Рис. 5. Распределения падающего радиационного теплового потока вдоль свободной поверхности (а), а также вдоль кварцевого (Ь) и графитового (с) тиглей Вертикальными пунктирами выделены особенности кривых, соответствующих конструкции тиглей (см. рис. 4)
На входной границе задается постоянное значение скорости прокачки газа V = 0,66 м/с. Скорости вращения тигля юс и кристалла ю^ составляли —5 и 20 об/мин соответственно. Свойства материалов, использованные в расчетах, представлены в табл. 2.
Определяющими параметрами данной задачи являлись число Прандтля Рг, число Грасгофа Ог, число Рэлея Ка, вращательные числа Рейнольдса Кес (тигля) и Ке^ (кристалла), число Марангони Ма, а также число ВК, характеризующее влияние сдвигового напряжения газа на течение расплава вдоль свободной поверхности. Они определяются следующим образом:
Таблица 2
Характеристики веществ, использованные в расчетах сопряженной задачи
Вещество Значение параметра
р, кг/м3 X, Вт/(м-К) с , Дж/(кгК) £ ,, гаа
Твердый кремний 2330 44-0,0138-Т 687-0,236-7 0,9016-0,00026208-7
Кварц 2650 4 1232 0,85
Графит 2000 70,7-0,0191-7 2019 0,80
Аргон Р0М V 0,01-2,5-10-57 532 -
И = 8,466-10-6 + 5,365-10-6-7 — 8,68240-127 Пах; р0 = 0,01 кг/м3; р0 = 3000 Па; М = 40 а.е.м.
Примечание. Свойства жидкого кремния уже приведены в табл. 1.
Обозначения: р0 — давление газа в ростовой установке; М - молекулярная масса; остальные соответствуют приведенным в табл. 1.
Рг = = 1,1 -10-2, (12)
X
г3
0г = ^—ь^к-= 3,4-108, (13)
V2
Ra = ОгРг = 3,7106, (14)
Я 2
Rec = = 4,9 .¡о4, (15)
V Я2
Re, = = 1,7 -104, (16)
V
Ма ^^^ = 6,7 -103, а7)
дТ Ьца
DN = \Т ^ ЯР= 2,6 -106. (18)
Ц
Здесь ЛТЫШ = 25,2 К — перепад температуры между точкой пересечения фронта кристаллизации с осью симметрии и точкой пересечения границы раздела расплав — тигель с осью симметрии; V, м2/с, — кинематический коэффициент вязкости, V = ц/р; Л 7^ = 19,3 К
— перепад температуры между тройными точками расплав-газ-кристалл и расплав-газ-тигель; Ь, м, — характерная протяженность свободной поверхности, Ь = Яс - Я^; а, м2/с,
— коэффициент температуропроводности, а = v/Pг; , Па, — среднее значение газового сдвигового напряжения вдоль свободной поверхности.
Для исследования сеточной чувствительности были проведены расчеты на трех сетках (с использованием модели Вольфштайна). Базовая сетка содержала около 20 тыс. ячеек, размер первой приповерхностной ячейки составлял около 0,3 мм, размер ячейки в объеме расплава — около 3,6 мм (см. рис. 4). Огрубленная и измельченная сетки получались путем уменьшения и увеличения числа ячеек в 2 раза по каждому направлению в расплаве и в 1 — 2 раза в остальных блоках.
Результаты расчетов сопряженной задачи
На рис. 6 приведены радиальные профили температуры, модуля меридиональной составляющей скорости на глубине 2 см от свободной поверхности расплава, а также распределения температуры вдоль границы раздела расплав — тигель, полученные в расчетах на трех сетках. Как видно из данных распределений, различия между результатами, полученными на базовой и измельченной сетках, меньше, чем при сравнении данных, полученных на огрубленной и базовой сетках, из чего можно сделать вывод о сеточной сходимости проведенных вычислений. Ввиду небольшой разницы между решениями, полученными на базовой и измельченной сетках, последующие расчеты были проведены на базовой сетке.
Рис. 6. Сравнение радиальных профилей температуры (а) и модуля меридиональной составляющей скорости (Ь) на глубине 2 см от свободной поверхности расплава с угловым распределением температуры вдоль границы раздела расплав — тигель (с). Также сравниваются
результаты использования огрубленной (1), базовой (2) и измельченной (3) сеток. Углы 0° и 90° соответствуют точкам пересечения стенки тигля с осью симметрии и со свободной
поверхностью расплава, соответственно
Рис. 7. Распределения температуры (верхний ряд), модуля меридиональной составляющей скорости (нижний ряд) по расплаву, а также векторные поля скорости (обозначены стрелками), полученные путем ILES- и RANS-расчетов по моделям Вольфштайна, Чена и Ментера
(слева направо)
На рис. 7 сравниваются поля температуры, модуля меридиональной составляющей скорости, а также векторов скорости, полученные с использованием ILES- и RANS-под-ходов. Перепад температуры, наиболее близкий к данным ILES-расчета, дает модель Вольфштайна, в то время как модель Чена предсказывает заниженный перепад, а модель Ментера, напротив, завышенный. Вихревая структура в окрестности вертикальной стенки тигля, отчетливо проявляющаяся в ILES-решении, воспроизводится и по моделям Вольфштайна и Ментера, тогда как модель Чена предсказывает иную структуру течения: с относительно низкими скоростями на периферии расплава.
Следует отметить, что все рассматриваемые в данной работе модели турбулентности предсказывают выраженное нисходящее течение расплава в окрестности оси симметрии, что противоречит результатам ILES-вычислений, где высоких скоростей расплава в под-кристальной области не наблюдается. Наиболее интенсивное нисходящее течение предсказывает модель Ментера, при этом она дает наименьший уровень турбулентной вязкости. _
Сравнение сдвиговой компоненты u'v' тензора рейнольдсовых напряжений, полученных в рамках вычислений ILES и RANS, представлено на рис. 8 (верхний ряд графиков). Здесь и' соответствует радиальной пульсации скорости, а_V,— осевой. Можно отметить качественное несоответствие распределений компоненты u'v', предсказываемых ILES- и RANS-подходами. В ILES-расчете наиболее высокие абсолютные значения u'v' наблюдаются в окрестности мениска, вблизи тройной точки расплав-газ-тигель и у закругленной части стенки тигля. Как установлено в работе [19], эта особенность связана с анизотропией пульсаций скорости, которая заключается в более сильном демпфировании нормальных пульсаций, по сравнению с продольными пульсациями у твердой стенки, а также демпфированием только нормальных пульсаций у свободной поверхности._
Таким образом, качественное различие распределений компоненты u'v' в ILES- и RANS-расчетах обусловлено отсутствием в концепции изотропной вязкости фактора, вызывающего поверхностную анизотропию.
В среднем и нижнем рядах графиков рис. 8 дается сравнение компонент вектора турбулентного теплового потока. ILES-подход предсказывает падение компоненты v'T' у горизонтальной части свободной поверхности, в то время как для компоненты и T' наблюдаются высокие отрицательные значения. Средний уровень турбулентного теплового потока согласуется в ILES- и RANS-расчетах. Ввиду отсутствия специальных поправок, учитывающих влияние свободной поверхности на турбулентность, уровень турбулентной вязкости на свободной поверхности получается того же порядка, что и в объеме расплава. Это, в свою очередь, приводит в RANS-расчетах к высоким значениям как горизонтальной, так и вертикальной компонент потока, что противоречит результатам ILES-расчета, где имеет
Рис. 8. Распределения сдвиговой компоненты u'v' тензора рейнольдсовых напряжений (верхний ряд) и компонент вектора турбулентного теплового потока (средний и нижний ряды), полученные в ILES- и RANS-расчетах по моделям Вольфштайна, Чена и Ментера
(слева направо)
место выраженная анизотропия теплопереноса, а снижение значений компоненты v'T' обусловлено демпфированием нормальных пульсаций скорости.
Заключение
Проведено сопоставление результатов RANS-расчетов турбулентной конвекции в расплаве кремния, полученных по нескольким моделям турбулентности с изотропной вязкостью, с ранее опубликованными данными ILES-вычислений для аналогичных условий.
Модели изотропной турбулентной вязкости в среднем отражают характеристики тепло-переноса в расплаве установок выращивания кристаллов по методу Чохральского, если основываться на температурных распределениях. Однако детали тепло- и массопереноса в расплаве, важные для изучения влияния ростовых параметров с целью оптимизации процесса выращивания, не могут быть корректно учтены в рамках моделей изотропной вязкости. Вместе с тем, можно сделать вывод об отсутствии преимуществ моделей с двумя дифференциальными уравнениями для предсказания структуры течения и поля температуры, по сравнению с моделью, включающей только одно уравнение. Среди трех рассмотренных RANS-моделей наиболее сильные различия в структуре конвекции, по сравнению с данными ILES, наблюдаются в случае k-ю SST-модели.
Для преодоления недостатков RANS-моделей с изотропной вязкостью можно использовать модели переноса рейнольдсовых напряжений, требующие дополнительного решения нескольких, существенно нелинейных дифференциальных уравнений. Однако это может приводить к численным сложностям для получения стационарного решения [24]. Другой путь — это проблемно-ориентированная модификация исходных моделей изотропной вязкости с алгебраическим введением факторов, которые, как показано в нашей работе [25], могут продуцировать желаемую анизотропию тензора рейнольдсовых напряжений и вектора турбулентного теплового потока, входящих в осредненные по Рейнольдсу уравнения движения и энергии. Представленное в настоящей статье сопоставление результатов RANS- и ILES-расчетов позволяет заключить, что, применительно к задачам расчета конвекции в тиглях установок метода Чохральского, в качестве исходной RANS-модели для модификации целесообразно взять либо однопараметрическую k-модель, либо двухпара-метрическую k-8 модель.
СПИСОК ЛИТЕРАТУРЫ
1. Vegad M., Bhatt N. M. Review of some aspects of single crystal growth using Czochralski crystal growth technique // Procedia Technology. 2014. Vol. 14. 2nd International Conference on Innovations in Automation and Mechatronics Engineering (ICIAME 2014). Pp. 438—446.
2. Wagner C., Friedrich R. Direct numerical simulation of momentum and heat transport in idealized Czochralski crystal growth configurations // International Journal of Heat and Fluid Flow. 2004. Vol. 25. No. 3. Pp. 431-443.
3. Raufeisen A., Breuer M., Botsch T., Delgado A. DNS of rotating buoyancy- and surface tension-driven flow // International Journal of Heat and Mass Transfer. 2008. Vol. 51. No. 25-26. Pp. 6219-6234.
4. Gräbner O., Mühe A., Müller G., Tomzig E., Virbulis J., Ammon W. V. Analysis of turbulent flow in silicon melts by optical temperature measurement // Materials Science and Engineering. B. 2000. Vol. 73. No. 1-3 Pp. 130-133.
5. Kobayashi S., Miyahara S., Fujiwara T., Kubo T., Fujiwara H. Turbulent heat transfer through the melt in silicon Czochralski growth // Journal of Crystal Growth. 1991. Vol. 109. No. 1-4. Pp. 149-154.
6. Kinney T. A., Brown R. A. Application of turbulence modeling to the integrated hydrodynamic thermal-capillary model of Czochralski crystal growth of silicon // Journal of Crystal Growth. 1993. Vol. 132. No. 3-4. Pp. 551-574.
7. Lipchin A., Brown R. A. Comparison of three turbulence models for simulation of melt convection in Czochralski crystal growth of silicon // Journal of Crystal Growth. 1999. Vol. 205. No. 1-2. Pp. 71-91.
8. Jones W. P., Launder B. E. The prediction of laminarization with a two-equation model of turbulence // International Journal of Heat and Mass Transfer. 1972. Vol. 15. No. 2. Pp. 301-314.
9. Verma S., Dewan A. Thermofluid characteristics of Czochralski melt convection using 3D URANS computations // ASME Journal of Thermal Science and Engineering Applications. 2019. Vol. 11. No. 6. P. 061017.
10. Launder B. E., Sharma B. I. Application of the energy-dissipation model of turbulence to the calculation of flow near a spinning disc // Letters in Heat and Mass Transfer. 1974. Vol. 1. No. 2. Pp. 131-137.
11. Menter F. R. Zonal two equation k-ю turbulence models for aerodynamic flows // Proceedings of the 24th Fluid Dynamics Conference. July 6 — 9, 1993, Orlando, Florida, USA. American Institute of Aeronautics & Astronautics (AIAA). 1993. Paper 1993—2906.
12. Menter F. R., Kuntz M., Langtry R. Ten years of industrial experience with the SST turbulence model // "Turbulence, Heat and Mass Transfer 4". Proceedings of the Symposium on Turbulence, Heat and Mass Transfer. 12 — 17 October, 2003. Antalya, Turkey. Ed. by Hanjalic K., Nagano Y., Tummers M. J. Book Series: Turbulence, Heat and Mass Transfer. Vol. 4. USA: Begell House Inc.,
2003. Pp. 625—632.
13. Lipchin A., Brown R. A. Hybrid finite-volume/finite-element simulation of heat transfer and melt turbulence in Czochralski crystal growth of silicon // Journal of Crystal Growth. 2000. Vol. 216. No. 1—4. Pp. 192—203.
14. Kalaev V. V., Lukanin D. P., Zabelin V. A., Makarov Y. N., Virbulis J., Dornberger E., Ammon W. V. Prediction of bulk defects in CZ Si crystals using 3D unsteady calculations of melt convection // Materials Science in Semiconductor Processing. 2002. Vol. 5. No. 4—5. Pp. 369—373.
15. Kalaev V. V., Lukanin D. P., Zabelin V. A., Makarov Y. N., Virbulis J., Dornberger E., Ammon W. V. Calculation of bulk defects in CZ Si growth impact of melt turbulent fluctuations // Journal of Crystal Growth. 2003. Vol. 250. No. 1—2. Pp. 203—208.
16. Wetzel T., Virbulis J., Muiznieks A., Ammon W. V., Tomzig E., Raming G., Weber M. Prediction of the growth interface shape in industrial 300 mm CZ Si crystal growth // Journal of Crystal Growth.
2004. Vol. 266. No. 1—3. Pp. 34—39.
17. Gräbner O., Müller G., Virbulis J., Tomzig E., Ammon W. V. Effects of various magnetic field configurations on temperature distributions in Czochralski silicon melts // Microelectronic Engineering. 2001. Vol. 56. No. 1—2. Pp. 83—88.
18. Müller G, Mühe A., Backofen R., Tomzig E., Ammon W. V. Study of oxygen transport in Czochralski growth of silicon // Microelectronic Engineering. 1999. Vol. 45. No. 2—3. Pp. 135—147.
19. Borisov D. V., Kalaev V. V. ILES of melt turbulent convection with conjugated heat transfer in the crucible and gas flow for Czochralski silicon crystal growth system // Journal of Crystal Growth. 2021. Vol. 573. 1 November. P. 126305.
20. Программный пакет CGSim [Электронный ресурс]. URL: http://www.str-soft.com/ products/CGSim/ (Дата обращения: 29.03.2022).
21. Wolfshtein M. The velocity and temperature distribution in one-dimensional flow with turbulence augmentation and pressure gradient // International Journal of Heat and Mass Transfer. 1969. Vol. 12. No. 3. Pp. 301—318.
22. Chien K.-Y. Predictions of channel and boundary-layer flows with a low-Reynolds number turbulence model // American Institute of Aeronautics & Astronautics (AIAA) Journal. 1982. Vol. 20. No. 1. Pp. 33—38.
23. Kalaev V., Sattler A., Kadinski L. Crystal twisting in Cz Si growth // Journal of Crystal Growth. 2015. Vol. 413. 1 March. Pp. 12—16.
24. Ristorcelli J. R., Lumley J. L. A second-order turbulence simulation of the Czochralski crystal growth melt: the buoyantly driven flow // Journal of Crystal Growth.1993. Vol. 129. No. 1—2. Pp. 249—265.
25. Kalaev V., Borisov D., Smirnov A. A modified hypothesis of Reynolds stress tensor modeling for mixed turbulent convection in crystal growth // Journal of Crystal Growth. 2022. Vol. 580. 15 February. P. 126464.
REFERENCES
1. Vegad M., Bhatt N. M., Review of some aspects of single crystal growth using Czochralski crystal growth technique, Procedia Technol. 14 (2nd Int. Conf. on Innovations in Automation and Mechatronics Engin. (ICIAME 2014)) (2014) 438—446.
2. Wagner C., Friedrich R., Direct numerical simulation of momentum and heat transport in idealized Czochralski crystal growth configurations, Int. J. Heat & Fluid Flow. 25 (3) (2004) 431—443.
3. Raufeisen A., Breuer M., Botsch T., Delgado A., DNS of rotating buoyancy- and surface tension-driven flow, Int. J. Heat & Mass Transfer. 51 (25-26) (2008) 6219-6234.
4. Gräbner O., Mühe A., Müller G., et al., Analysis of turbulent flow in silicon melts by optical temperature measurement, Mater. Sci. Eng. B. 73 (1-3) (2000) 130-133.
5. Kobayashi S., Miyahara S., Fujiwara T., et al., Turbulent heat transfer through the melt in silicon Czochralski growth, J. Cryst. Growth. 109 (1-4) (1991) 149-154.
6. Kinney T. A., Brown R. A., Application of turbulence modeling to the integrated hydrodynamic thermal-capillary model of Czochralski crystal growth of silicon, J. Cryst. Growth. 132 (3-4) (1993) 551-574.
7. Lipchin A., Brown R. A., Comparison of three turbulence models for simulation of melt convection in Czochralski crystal growth of silicon, J. Cryst. Growth. 205 (1-2) (1999) 71-91.
8. Jones W. P., Launder B. E., The prediction of laminarization with a two-equation model of turbulence, Int. J. Heat & Mass Transfer. 15 (2) (1972) 301-314.
9. Verma S., Dewan A., Thermofluid characteristics of Czochralski melt convection using 3D URANS computations, ASME J. Therm. Sci. Eng. Appl. 11 (6) (2019) 061017.
10. Launder B. E., Sharma B. I., Application of the energy-dissipation model of turbulence to the calculation of flow near a spinning disc, Lett. Heat & Mass Transfer. 1 (2) (1974) 131-137.
11. Menter F. R., Zonal two equation k-© turbulence models for aerodynamic flows, Proc. 24th Fluid Dynamics Conf., July 6 - 9, 1993, Orlando, Florida, USA, American Institute of Aeronautics & Astronautics (AIAA), 1993. Paper 1993-2906.
12. Menter F. R., Kuntz M., Langtry R., Ten years of industrial experience with the SST turbulence model, In book: "Turbulence, Heat and Mass Transfer 4". Proc. Symp. on Turbulence, Heat & Mass Transfer. 12-17 Oct., 2003, Antalya, Turkey. Ed. by Hanjalic K., Nagano Y., Tummers M. J., Book Ser. Turbulence, Heat & Mass Transfer. Vol. 4, Begell House Inc., USA (2003) 625-632.
13. Lipchin A., Brown R. A., Hybrid finite-volume/finite-element simulation of heat transfer and melt turbulence in Czochralski crystal growth of silicon, J. Cryst. Growth. 216 (1-4) (2000) 192-203.
14. Kalaev V. V., Lukanin D. P., Zabelin V. A., et al., Prediction of bulk defects in CZ Si crystals using 3D unsteady calculations of melt convection, Mater. Sci. Semicond. Proc. 5 (4-5) (2002) 369-373.
15. Kalaev V. V., Lukanin D. P, Zabelin V. A., et al., Calculation of bulk defects in CZ Si growth impact of melt turbulent fluctuations, J. Cryst. Growth. 250 (1-2) (2003) 203-208.
16. Wetzel T., Virbulis J., Muiznieks A., et al., Prediction of the growth interface shape in industrial 300 mm CZ Si crystal growth, J. Cryst. Growth. 266 (1-3) (2004) 34-39.
17. Gräbner O., Müller G., Virbulis J., et al., Effects of various magnetic field configurations on temperature distributions in Czochralski silicon melts, Microelectr. Eng. 56 (1-2) (2001) 83-88.
18. Müller G, Mühe A., Backofen R., et al., Study of oxygen transport in Czochralski growth of silicon, Microelectr. Eng. 45 (2-3) (1999) 135-147.
19. Borisov D. V., Kalaev V. V., ILES of melt turbulent convection with conjugated heat transfer in the crucible and gas flow for Czochralski silicon crystal growth system, J. Cryst. Growth. 573 (1 November) (2021) 126305.
20. Software package CGSim [Electronic resource]. URL: http://www.str-soft.com/products/ CGSim/ (Last accessed date: March 29, 2022).
21. Wolfshtein M., The velocity and temperature distribution in one-dimensional flow with turbulence augmentation and pressure gradient, Int. J. Heat & Mass Transfer. 12 (3) (1969) 301-318.
22. Chien K.-Y., Predictions of channel and boundary-layer flows with a low-Reynolds number turbulence model, AIAA J. 20 (1) (1982) 33-38.
23. Kalaev V., Sattler A., Kadinski L., Crystal twisting in Cz Si growth, J. Cryst. Growth. 413 (1 March) (2015) 12-16.
24. Ristorcelli J. R., Lumley J. L., A second-order turbulence simulation of the Czochralski crystal growth melt: the buoyantly driven flow, J. Cryst. Growth. 129 (1-2) (1993) 249-265.
25. Kalaev V., Borisov D., Smirnov A., A modified hypothesis of Reynolds stress tensor modeling for mixed turbulent convection in crystal growth, J. Cryst. Growth. 580 (15 February) (2022) 126464.
СВЕДЕНИЯ ОБ АВТОРАХ
БОРИСОВ Дмитрий Витальевич — инженер-программист АО «Группа СТР» — ООО «Со-фт-Импакт», Санкт-Петербург, Россия.
194044, Россия, г. Санкт-Петербург, Большой Сампсониевский пр., 64
dmitriy.borisov@str-soft.com
ORCID: 0000-0002-7481-8265
КАЛАЕВ Владимир Владимирович — кандидат физико-математических наук, технический директор АО «Группа СТР» — ООО «Софт-Импакт», Санкт-Петербург, Россия. 194044, Россия, г. Санкт-Петербург, Большой Сампсониевский пр., 64 vladimir.kalaev@str-soft.com ORCID: 0000-0001-9231-0740
THE AUTHORS
BORISOV Dmitry V.
STR Group, Inc. — Soft-Impact, Ltd.
64, Bolshoi Sampsonievskii Ave., St. Petersburg, 194044, Russia
dmitriy.borisov@str-soft.com
ORCID: 0000-0002-7481-8265
KALAEV Vladimir V.
STR Group, Inc. - Soft-Impact, Ltd.
64, Bolshoi Sampsonievskii Ave., St. Petersburg, 194044, Russia
vladimir.kalaev@str-soft.com
ORCID: 0000-0001-9231-0740
Статья поступила в редакцию 29.03.2022. Одобрена после рецензирования 30.05.2022. Принята 30.05.2022.
Received 29.03.2022. Approved after reviewing 30.05.2022. Accepted 30.05.2022.
© Санкт-Петербургский политехнический университет Петра Великого, 2022
