Содержание к диссертации
Введение
1. Методы численного моделирования турбулентных термоконвективных течений: обзор 12
1.1. Предварительные замечания 12
1.2. Метод прямого численного моделирования турбулентности 13
1.3. Решение осредненных по Рейнольдсу уравнений переноса 18
1.3.1. Уравнения Рейнольдсаи подходы к их замыканию 18
1.3.2. Модели турбулентной вязкости 20
1.3.2.1. Исходные понятия и классификация 20
1.3.2.2. Алгебраические модели 21
1.3.2.3. Модели с одним уравнением 23
1.3.2.4. Модели с двумя уравнениями 26
1.3.3. Заметки об эволюции численных подходов к описанию турбулентной конвекции Релея-Бенара в рамках RANS 27
1.4. Метод моделирования крупных вихрей 30
1.4.1. Операция пространственной фильтрации 31
1.4.2. Пространственно отфильтрованные определяющие уравнения 34
1.4.3. Модели подсеточных напряжений 35
1.4.3.1. Алгебраические модели типа модели Смагоринского 35
1.4.3.2. Динамические модели 39
1.4.3.3. Модели близких масштабов и смешанные модели 42
1.4.3.4. Модели с одним уравнением 43
1.4.3.5. Априорное тестирование SGS моделей 44
1.4.4. Моделирование пристенных областей в методе LES 44
1.4.5. Заключительные замечания относительно метода LES 46
1.5. Комбинированные подходы к моделированию турбулентности 47
1.5.1. Зональный метод 48
1.5.2. Метод встраивания зоны RANS в пристенную часть сетки LES 49
1.5.3. Модели, основанные на сравнении или взвешивании турбулентной вязкости 50
1.5.4. Метод моделирования отсоединенных вихрей 51
2. Математическая модель и численный метод 55
2.1. Математическая модель 55
2.1.1. Определяющие уравнения 55
2.1.2. Пристенные функции 57
2.2. Численный метод 59
2.2.1. Предварительные замечания 59
2.2.2. Запись определяющих уравнений в обобщенной системе координат 60
2.2.3. Численная схема 61
2.2.4. Предшествующая верификация программного комплекса 64
2.3. Реализация замыкающих моделей LES и RANS/LES в программном комплексе SINF 66
2.3.1. Реализация моделей LES 66
2.3.2. Модели RANS/LES и их реализация 68
3. Турбулентная конвекция Релея-Бенара в областях простой формы 72
3.1 Обзор аналитических и экспериментальных работ 72
3.1.1. Возникновение конвекции. Переходные режимы 72
3.1.2. Аналитические подходы к описанию турбулентной конвекции Р-Б 73
3.1.3. Экспериментальное изучение турбулентной конвекции Р-Б 76
3.2. Турбулентная конвекция во вращающемся горизонтальном слое, подогреваемом снизу 80
3.2.1. Предварительные замечания 80
3.2.2. Постановка и вычислительные аспекты задачи. Варианты расчетов 83
3.2.3. Структура турбулентной конвекции во вращающемся и неподвижном слое 87
3.2.3.1. Эффекты вращения и влияние числа Релея 87
3.2.3.2. Сопоставление результатов, полученных для Ra = 1.13x10s
при разных моделях турбулентности 94
3.2.4. Статистические и спектральные характеристики. Масштабные закономерности 94
3.2.5. Предсказание теплопроводящих свойств слоя: числаNu 105
3.2.6. Заключительные замечания 108
3.3. Турбулентная конвекция ртути в подогреваемой снизу цилиндрической полости 108
3.3.1. Предварительные замечания 108
3.3.2. Постановка и вычислительные аспекты задачи 115
3.3.3. Общая структура течения. Влияние числа Релея. Глобальная циркуляция в полости 118
3.3.4. Характеристики осредненных полей. Толщины пограничных слоев 129
3.3.5. Спектральные характеристики конвекции Анализ турбулентных пульсаций 133
3.3.5. Интегральные тепловые потоки. Зависимость Nil (Ra) 137
3.4. Турбулентная конвекция воды в подогреваемой снизу кубической полости... 140
3.4.1. Предварительные замечания 140
3.4.2. Постановка и вычислительные аспекты задачи 146
3.4.3. Структура Р-Б конвекции в кубической полости 148
3.4.3. Анализ турбулентных пульсаций. Интегральные тепловые потоки 154
4. Турбулентная конвекция расплава кремния в емкостях с геометрией, типичной для тиглей метода Чохральского 158
4.1. Введение в проблему 158
4.2. Постановка и вычислительные аспекты задачи 162
4.3. Прямое численное моделирование конвекции расплава и транспорта кислорода при Ra = 5х105 166
4.3.1. Предварительные замечания 166
4.3.2. Влияние вращения тигля на конвекцию расплава и концентрацию кислорода 167
4.4. Моделирование конвекции расплава и транспорта кислорода при Ra = 8.2x10 175
4.4.1. Результаты расчетов конвекции расплава методом LES 175
4.4.1.1. Влияние скорости вращения тигля на конвекцию расплава 176
4.4.1.2. Влияние расхода аргона на конвекцию расплава 183
4.4.2. Результаты расчетов конвекции расплава методом RANS/LES 185
Заключение 190
Литература 192
Приложение Ш
- Метод прямого численного моделирования турбулентности
- Комбинированные подходы к моделированию турбулентности
- Реализация замыкающих моделей LES и RANS/LES в программном комплексе SINF
- Общая структура течения. Влияние числа Релея. Глобальная циркуляция в полости
Введение к работе
Изучение конвективных движений жидкости, возникающих вследствие пространственной неоднородности плотности в полях массовых сил, имеет глубокое фундаментальное и важное практическое значение для самых разнообразных областей науки и техники таких, как геофизика, астрофизика, теплоэнергетика, микроэлектроника, кондиционирование и вентиляция, химические технологии. Неудивительно, что разнообразные формы конвективных течений постоянно привлекают повышенный интерес исследователей. Достижения в этой области отражены в ряде монографий и огромном числе отдельных научных публикаций [Гершуни, 1973; Джалурия, 1983; Роди, 1984; Монин, 1988; Сб. Гидромеханика..., 1990; Гебхарт, 1991]. Вместе с тем и по сегодняшний день актуальной остается проблема расчетного предсказания характеристик турбулентных режимов конвекции, представляющих наибольший практический интерес.
В ряду высокотехнологичных практических областей, в которых особенно востребованы углубленные знания о турбулентной конвекции, важное место занимает процесс выращивания кристаллов полупроводников из расплава по методу Чохральского [Полежаев, 1984; Мюллер, 1991]. Имеющийся на сегодняшний день объем сведений, в том числе фундаментального характера, о свойствах режимов течения полупроводниковых расплавов крайне недостаточен. Проведение экспериментов, направленных на изучение турбулентных режимов конвекции в приближенных к практике условиях является весьма трудоемкой, сопряженной со значительными материальными и временными затратами, а зачастую и нереализуемой задачей. В связи с этим численное моделирование турбулентной конвекции рассматривается сегодня как реальная альтернатива экспериментальному подходу.
В определенной степени настоящая диссертационная работа направлена на создание адекватных методов расчета турбулентных режимов конвекции расплава в установках метода Чохральского, в частности, расплава кремния. Очевидно, однако, что при отработке методов численного моделирования сложных конвективных движений, имеющих место в реальных практических установках, целесообразно опираться на опыт, приобретенный в работе с модельными задачами, предполагающими использование упрощенной геометрии и идеализированных граничных условий. Это открывает возможность верификации расчетов путем сопоставления с данными, полученными в результате проведения сравнительно недорогих лабораторных экспериментов. Рассматриваемое перед обращением к проблемам технологии Чохральского свободно-конвективное движение жидкости в подогреваемых снизу областях (конвекция Релея Бенара), является удобной моделью для отработки алгоритмов и методов численного описания структуры турбулентной конвекции, а также процессов тепло и массопереноса. Кроме того, данная задача сама по себе имеет широкое прикладное значение, и, обладая сложной внутренней структурой, представляет интерес с общих позиций [Буссе, 1983; Зимин, 1988; Siggia, 1994].
Несмотря на колоссальный прогресс в вычислительных ресурсах, а также достигнутые в последнее время значительные успехи в области построения эффективных численных алгоритмов решения задач гидродинамики и тепломассообмена, проблемы, возникающие при численном моделировании турбулентных конвективных течений по-прежнему остаются трудноразрешимыми. Среди трех используемых сегодня подходов к численному описанию турбулентной конвекции наибольшей привлекательностью, в силу своей строгости и надежности результатов расчетов при выполнении условий полной разрешимости всех составляющих движения, обладает метод прямого численного моделирования турбулентности {Direct Numerical Simulation, DNS). Однако возможности по применению DNS с полным разрешением ограничиваются модельными задачами в упрощенной геометрии при относительно невысоких числах Рейнольдса (Релея), что связано с очень быстрым ростом требований по вычислительным ресурсам при попытках продвинуться вверх по числам Рейнольдса. Метод решения осредненных по Рейнольдсу уравнений Навье-Стокса {Reynolds-Aver aged Navier-Stokes, RANS), требующий значительно меньших вычислительных затрат, чем DNS, до сих пор остается наиболее распространенным подходом к моделированию турбулентных течений. Вместе с тем, общепризнанно, что результаты расчетов по методу RANS очень чувствительны к выбору той или иной модели турбулентности, а иногда и просто не способны отразить характерные особенности, присущие термоконвективным течениям. Метод моделирования крупных вихрей {Large Eddy Simulation, LES) предполагает выполнение точного расчета переноса импульса и тепла лишь крупными, энергетически важными структурами и позволяет проводить расчеты при значительно более высоких числах Рейнольдса (Релея), в сравнении с DNS, с привлечением более простых замыкающих моделей, чем те, которые используются в подходах RANS [Eidson, 1985; Wong, 1993; Piomelli, 1998; Frohlich, 2000]. Особая привлекательность метода LES применительно к расчетам термоконвективных задач объясняется способностью адекватно воспроизводить эволюцию во времени вихревых структур, которые в значительной мере определяют свойства таких течений. Хорошо известно, что наибольшие трудности метода LES связаны с расчетом пристенных зон течения. Преодолеть эти трудности, проявляющиеся при моделировании пристенных областей потока, помогают развиваемые в последнее время гибридные под - 10 ходы, сочетающие в себе те или иные элементы методов RANS, LES и DNS [Balaras, 1996; Spalart, 2000; Smirnov, 20026].
Исходя из изложенных соображений, определены основные цели работы. Диссертационная работа направлена на
1) реализацию ряда замыкающих моделей метода моделирования крупных вихрей (LES) и гибридных подходов RANS/LES для численных расчетов турбулентных течений на основе базового программного комплекса SINF, развиваемого сотрудниками кафедры гидроаэродинамики СПбГПУ;
2) проведение методических расчетов, направленных на апробацию внедренных моделей в рамках методов LES и RANS/LES, на примере задачи о турбулентной конвекции в подогреваемом снизу вращающемся слое модельной жидкости (Pr = 1);
3) проведение численного моделирования на основе методов DNS и RANS/LES турбулентных режимов конвекции воды (Pr = 7) и ртути (Pr = 0.025) в подогреваемых снизу замкнутых полостях простой формы, в диапазоне чисел Релея 105 Ra 5х109;
4) численное исследование с использованием методов DNS, LES и RANS/LES турбулентной конвекции расплава кремния в емкости, геометрия которой типична для тиглей метода Чохральского, включая вопросы транспорта выделяющегося на стенках тигля кислорода.
Первая глава диссертации посвящена обзору методов численного моделирования турбулентных термоконвективных течений. В ней даются оценки требуемых вычислительных затрат при расчетах турбулентной конвекции по методу DNS. Приводится общепринятая классификация моделей метода RANS, анализируются проблемы, возникающие при расчетах на их основе задач турбулентной термоконвекции, в частности конвекции Релея-Бенара. Рассматриваются вопросы пространственной фильтрации в методе LES, приводится обзор моделей подсеточных масштабов. Излагаются подходы к моделированию пристенных областей в рамках метода LES, а также подходы, основанные на комбинировании методов RANS и LES.
Во второй главе дается изложение математической модели и численного метода, на основе которого были проведены расчеты (задачи турбулентной конвекции рассматриваются в работе в рамках модели несжимаемой ньютоновской жидкости с постоянными физическими свойствами, для учета эффектов плавучести в поле силы тяжести используется приближение Буссинеска). Излагаются вопросы реализации замыкающих моделей метода LES (модели Смагоринского и модели с одним уравнением) и гибридных моделей RANS/LES (алгебраической и модели с одним уравнением) в использованном для расчетов программном "конечно-объемном" комплексе SINF, второ -11 го порядка точности по времени и пространству. Применительно к методу LES рассматриваются способы демпфирования подсеточной вязкости у стенок, а также особенности реализации пристенных функций.
Третья глава посвящена численному моделированию турбулентной конвекции Ре-лея-Бенара, развивающейся в областях простой геометрической формы. Проводится дополнительная верификация программного комплекса SINF для расчетов турбулентной конвекции по методу DNS при умеренных числах Релея. С целью начальной отработки реализованных LES и RANS/LES моделей расчета турбулентных термоконвективных задач рассматривается течение модельной жидкости (7V=1) во вращающемся горизонтальном слое, подогреваемом снизу, при Rax 10 . Путем сопоставления с известными из литературы высокоточными данными DNS выявляются некоторые слабые и сильные стороны использованных моделей LES и RANS/LES применительно к задачам рассматриваемого класса, а также дается ряд рекомендаций по их дальнейшему использованию в расчетах. Представляются результаты численного моделирования на основе метода RANS/LES турбулентной конвекции в подогреваемой снизу кубической полости, доверху наполненной водой (Рг = 7), и конвекции ртути (Pr = 0.025) в цилиндрической полости, при варьировании числа Релея в диапазоне 108 і?а 5х109. Анализируются специфические черты турбулентной конвекции, такие как эволюция термиков, глобальная циркуляция жидкости в полости. Проводится сопоставление с известными данными экспериментов рассчитанных толщин пограничных слоев и чисел Нуссельта.
В четвертой главе рассматриваются вопросы численного моделирования турбулентной конвекции расплава кремния в емкости с геометрией, типичной для тиглей метода Чохральского при числе Релея Ra = 8.2х106, характеризующем реальные условия роста, и при уменьшенном числе Ra = 5x105. Представляются детальные данные о пульсирующих полях скорости, температуры и концентрации кислорода при различных значениях определяющих параметров. Анализируются эффекты, связанные с изменением скорости вращения тигля и условий на свободной поверхности. Для режима с Ra = 8.2x106 проводится сопоставление с доступными экспериментальными данными по средним полям и пульсациям температуры в расплаве.
В заключении представлены основные результаты и выводы, полученные в работе.
Метод прямого численного моделирования турбулентности
Метод прямого численного моделирования турбулентности (в иностранной литературе - Direct Numerical Simulation, DNS) является наиболее привлекательным и надежным подходом, предоставляя возможность всестороннего изучения турбулентных явлений в случаях, когда экспериментальные исследования затруднены или вообще невозможны [Grotzbach, 1999]. Этот подход является максимально строгим, так как базируется лишь на одном, достаточно обоснованном предположении о применимости уравнений Навье-Стокса для представления турбулентного движения среды. В приложении к термоконвективным течениям метод DNS предполагает решение дискретизированных уравнений (1.1)-(1.3) без привлечения каких-либо дополнительных физических моделей. Для получения аккуратного решения в рамках такого подхода все значимые временные и пространственные масштабы течения должны разрешаться полностью. Отношение наибольших и наименьших пространственных масштабов течения определяет число степеней свободы, необходимое для численного представления поля течения в любом из трех измерений. Наибольший масштаб турбулентных образований в ограниченной области имеет порядок характерного размера этой области L. В качестве наименьшего масштаба длины обычно выбирается колмогоровский масштаб [Колмогоров, 1941] где є - скорость диссипации энергии турбулентных пульсаций, оцениваемая как = V /L. Тогда для отношения наибольшего к наименьшему масштабу можно получить следующую оценку [Шуманн, 1984]: Как видно, величина этого отношения растет с увеличением числа Рейнольдса. Требуемое для полного разрешения в рамках DNS число расчетных узлов (контрольных объемов), играющих роль степеней свободы в конечно-разностных (конечно-объемных) схемах, с учетом трех пространственных направлений пропорционально Re9 4. Для получения статистически установившегося решения расчеты должны продолжаться в течение промежутка времени, по крайней мере, порядка L /V/,, а шаг по времени, в свою очередь, не должен превышать величину т] /Vb.
Таким образом, про " /А ведение полного расчета предполагает выполнение порядка Re шагов по времени. Следовательно, суммарная "стоимость" расчетов для "полного" DNS нарастает как Re3 или Ra3 2. При изучении термоконвективных задач к требованию по разрешению вихрей колмогоровских масштабов г] добавляется дополнительное условие точного численного представления так называемых диффузионных тепловых масштабов Бэтчелора [Kerr, 1996] Понятно, что при численном моделировании конвективных течений жидкостей, характеризующихся высокими значениями числа Прандтля, требования по разрешению в рамках DNS становятся еще более жесткими, поскольку минимальным масштабом в данном случае становится 7]ь. Для оценки качества результатов DNS в задачах о конвекции Релея-Бенара (Р-Б) часто пользуются соотношениями, которые определяют требующиеся для разрешения минимальных масштабов размеры ячеек расчетной сетки в центральной части области, где предполагается изотропность линейных размеров ячеек 5С« VCV3, Vc - объем ячейки. Согласно оценкам [Grotzbach, 1983], условия разрешения удаленных от стенок областей выполнены, если где Nu - число Нуссельта. К условиям (1.7) необходимо добавить требование качественного разрешения высокоградиентных слоев у горизонтальных стенок области течения. В ряде работ по численному моделированию турбулентной конвекции Р-Б [Grotzbach, 1983; Verzicco, 1999] приводятся оценки, в соответствии с которыми, для получения адекватных результатов DNS эти динамически важные области течения должны быть "покрыты", по меньшей мере, 5-7 узлами расчетной сетки. При этом максимальное удлинение ячеек у стенок не должно превышать 10. Предварительная приближенная оценка толщин температурного и скоростного пристенных слоев может быть осуществлена на основе соотношений:
Опираясь на соотношения (1.7), (1.8) можно указать, в качестве примера, требования по разрешению в задаче конвекции Р-Б в замкнутой полости с отношением ее высоты к ширине, равным единице. Если принять значения определяющих параметров Ra= 109, Рг=\, то в этом случае по данным экспериментов число Нуссельта приближенно равно 75. Тогда Sc 1.9x10 , Хт Лу 6.7x10" . Тогда минимально необходимое число расчетных узлов в вертикальном направлении получается приблизительно равным N «600; далее примем Nx = Ny = 2/3Nz. В итоге, требуемая размерность расчетной сетки составляет около 100 миллионов ячеек. Аналогичные вычисления для Ra = 108 дают требуемую размерность - 12 миллионов ячеек. Принимая во внимание произведенные оценки, можно констатировать, что на сегодняшний день возможности по применению DNS ограничиваются задачами при числах Релея, не превышающих 108 [Smirnov, 20026]. Тем не менее, высокоточные результаты прямого численного моделирования турбулентности даже при сравнительно невысоких значениях числа Релея являются весьма востребованными. Детальный анализ результатов таких расчетов может способствовать прояснению физических механизмов, определяющих специфику термоконвективных течений, послужить надежной опорой при разработке и тестировании подходов, позволяющих осуществить в вычислительном смысле "продвижение" вверх по числам Релея. Кроме того, некоторые классы термоконвективных течений характеризуются установлением "автомодельных" по числу Ra режимов течения уже при сравнительно невысоких значениях последнего, что дает возможность обобщения, и оправдывает экстраполяцию полученных результатов на режимы более развитой конвекции. Первые численные исследования переходных и слаботурбулентных режимов конвекции Релея-Бенара в трехмерной постановке относятся к началу-середине 1980-х годов [Grotzbach, 1982; Curry, 1984; Domaradzki, 1988]. Расчеты проводились
Комбинированные подходы к моделированию турбулентности
Наиболее прямой подход к комбинированию RANS и LES состоит в предварительном, до проведения расчетов разделении расчетной области на статические "зоны ответственности" RANS и LES. При выполнении расчетов в каждой из зон "работает" один из заранее назначенных для нее методов (RANS, либо LES). При этом указанные методы могут быть основаны на различных, не связанных между собой моделях турбулентности. Примеры использования зонального метода можно найти, в частности, в работах [Davidson, 2001; Temmerman, 2002; Hamba, 2002], где он привлекался для моделирования турбулентного течения в канале. В первой из упомянутых работ [Davidson, 2001] в зонах RANS использовалась "низкорейнольдсовая" к-со модель, а в зонах LES - модель с одним уравнением для ksgs [Yoshizawa, 1982]. В качестве характерного масштаба длины в зоне LES выбирался минимальный из линейных размеров ячейки расчетной сетки A = min(Ai), / = 1,2,3. В работе [Temmerman, 2002] метод RANS основывался на "низкорейнольдсовой" модели с одним уравнением для кинетической энергии турбулентности [Wolfstein, 1969], а метод LES - на алгебраической модели типа Смагоринского. Во всех расчетах в некоторой "буферной" области, включающей в себя точку переключения между зонами, наблюдалось почти ступенчатое изменение средней скорости (в несколько единиц и+) с переходом профиля скорости с прямой, которая соответствует логарифмическому закону (предсказываемому моделью RANS) на расположенную выше первой логарифмическую прямую, которую дает подсеточная модель.
При этом в области LES значения и превышали данные экспериментов, следствием чего явилось занижение в среднем на 15% значения коэффициента трения на стенке. Судя по литературе, отмеченная проблема почти всегда проявляется при попытках описания пристенных турбулентных течений на основе комбинирования подходов с радикально различающимся механизмом генерации турбулентных напряжений [Baggett, 1998; Piomelli, 2002]. Предположив, что основной причиной образования несогласованного профиля скорости является излишне резкое изменение вихревой вязкости в окрестности места переключения, что влечет за собой неаккуратное определение конвективных членов в уравнении баланса количества движения, Хамба [Hamba, 2002] предложил вычислительную технику, которая позволила уменьшить степень рассогласования на профиле и . В общих чертах техника основана на привлечении в расчетах двух наборов компонент скорости, соответствующих двум различным значениям ширины фильтра на границе переключения. В качестве моделей турбулентности были использованы низкорейнольдсовая к-є модель в области RANS и модель с одним уравнением [Yoshizawa, 1982] в области LES. Координата точки переключения у входила в параметры задачи, также как и ширина некоторой "буферной области", включающей значение у . Для обеспечения гладкого перехода между областями RANS и LES величина скорости диссипации турбулентной энергии в буферной области вычислялась как Следует заметить, что кроме априорных знаний о структуре течения, позволяющих заранее подходящим образом разбить расчетную область на статические RANS - LES - зоны, рассмотренный метод предполагает постановку специальных граничных условий для характеристик турбулентности (и/или вихревой вязкости) на границе раздела зон. 1.5.2. Метод встраивания зоны RANS в пристенную часть сетки LES
Данный метод, именуемый также "двухслойной моделью" (Two-Layer Model, TLM) турбулентности [Balaras, 1996; Cabot, 1996], нацелен на более аккуратное (по сравнению с методом пристенных функций) вычисление напряжения трения на стенке при расчетах по методу LES на сетках, не имеющих выраженного сгущения у стенок. При этом только на участке от стенки до ближайшего к ней узла (составляющем зону RANS), расчетная (базовая) LES - сетка дробится в вертикальном направлении. На этой RANS - сетке, "встроенной" в пристенную часть LES - сетки, численно решаются дифференциальные уравнения движения, замкнутые по какой-либо модели турбулентности RANS. В работе [Balaras, 1996] предлагается для сформированной RANS - зоны решать систему трехмерных нестационарных уравнений "тонкого пристенного слоя" (Thin Boundary Layer Equations, TELE) для параллельных стенке компонент скорости «f (/=1,3): и граничными условиями: u,(ym) = Umi, и/0,) = 0. Первое из граничных условий обеспечивает непрерывность параллельных стенке компонент скорости на границе раздела RANS и LES зон (при у = ут), сохраняя тем самым связь с актуальным решением во внешней LES - зоне. Турбулентная вязкость vt в зоне RANS вычисляется на основе характерных для RANS алгебраических соотношений с привлечением демпфирующих соотношений [Balaras, 1996] Величина напряжения трения на стенке, которая используется при решении LES - уравнений, определяется по градиентному соотношению, использующему най денное в результате решения (1.60) поле скорости: Для снижения вычислительных затрат в [Cabot, 1999] предлагается использовать упрощенный вид системы (1.60), в котором левая часть уравнений полагается равной нулю. В результате требуется решение двух простых, несвязанных друг с другом обыкновенных дифференциальных уравнений. Другой возможный подход предполагает использование аппроксимации некоторых членов из левой части уравнений (1.60) значениями величин во внешнем потоке [Hoffmann, 1995]. турбулентной вязкости Принципиально иной способ комбинирования RANS и LES предполагает использование единой сетки. Во всей расчетной области на основе какой-либо модели метода RANS и модели метода LES вычисляются два значения турбулентной вязкости (у т и vEES). После чего для определения величины турбулентной вязкости ц, входящей в уравнения движения, применяются два способа. В первом случае в качестве vt принимается минимальное из значений yNS и vEES. Такой подход был, в частности, использован в работе [Schumann, 1975]. Рассматривая задачу о турбулентном течении в канале, автор предложил разложить тензор напряжений ц- на "изотропную" часть (для представления которой использовалась SGS модель Смагоринского) и "анизотропную" часть, отражающую глобальный градиент скорости в потоке. "Анизотропная" часть тензора напряжений моделируется с использованием следующего алгебраического соотношения для вихревой вязкости:
Реализация замыкающих моделей LES и RANS/LES в программном комплексе SINF
Представленная в разделе 2.2 базовая версия ПК SINF была дополнена в ходе настоящей работы возможностью расчета турбулентных течений на основе методов LES и RANS/LES. В данном разделе излагаются отдельные аспекты реализации внедренных в ПК замыкающих моделей, а также пристенных функций, использованных в расчетах по методу LES. В рамках метода конечных объемов, на котором основана дискретизация пространственных операторов уравнений сохранения в ПК SINF, операция пространственной фильтрации метода LES может быть интерпретирована как осреднение по объемам ячеек расчетной сетки. В качестве локального масштаба длины в этом случае традиционно (также и в настоящей работе) используется кубический корень из объема ячейки расчетной сетки (см. (1.45)). Отметим, что требующиеся в моделях подсеточной турбулентности значения компонент тензора скоростей деформации в ПК SINF находятся в центрах ячеек путем применения интегральной формулировки для вычисления градиента скалярных функций (в данном случае - декартовых компонент скорости), при этом используется линейная интерполяция на грани ячейки компонент скорости из соседних ячеек. Применительно к методу LES никакой дополнительной пространственной фильтрации поля скорости с целью получения модуля тензора скоростей деформации входящего в алгебраическую модель Смагоринского (1.37), не производилось. Как следствие, реализация этой алгебраической модели подсеточной турбулентности оказалась относительно простой задачей. Математическая формулировка реализованной модели с одним уравнением, которая основана на решении уравнения для кинетической энергии неразрешаемого движения, выглядит следующим образом: Постоянным модели "по умолчанию" присваиваются следующие значения: С = 0.12, Сє= 0.7, ак - 1. Внедрение модели в ПК SINF производилось по примеру реализованной в нем высокорейнольдсовой к-є модели турбулентности RANS.
Для уменьшения подсеточной вязкости вблизи стенок может использоваться (при выборе соответствующей опции) демпфирующий множитель вида [Зябриков, 1987] при этом найденная на основе либо модели Смагоринского, либо модели одного уравнения величина vt умножается на D{y+). Привлечение демпфирующих соотношений подразумевает вычисление для всех внутренних узлов расчетной сетки расстояния до ближайшей стенки. В реализованном на сегодняшний день в ПК SINF алгоритме это расстояние находится путем отчасти оптимизированной процедуры перебора расстояний от данной внутренней точки до всех узлов, расположенных на стенках. На сетках большой размерности реализованный алгоритм требует существенных временных затрат (например, в расчетах Р-Б конвекции из раздела 3.2 настоящей работы, проведенных на сетке в 180,000 ячеек, при использовании процессора Intel Pentium-Ill, 1000 МГц, процесс вычисления расстояний до стенки занял 14 минут). Однако поскольку рассматриваемая операция является разовой, выполняемой до проведения основных (нестационарных) расчетов, требующих нескольких дней или недель счета, то в целом затраты на расчеты наименьшего расстояния до стенки оказываются непринципиальными. Для определения параметров течения вблизи стенки в рамках техники пристенных функций в ПК SINF реализован следующий алгоритм. На каждой итерации по времени установления для всех ближайших к стенке ячеек величина отношения uTIU вычисляется непосредственно по формулам (2.7). При этом ву+ из правой части (2.7) входит величина динамической скорости с предыдущей итерации. Дополнительно используется релаксационное соотношение где индекс М-1 соответствует величинам, полученным на предыдущей итерации. Величина трения на стенке определяется по формуле Релаксационное соотношение типа (2.37) используется и при вычислении компонент напряжения трения на стенке. В проведенных в настоящей работе расчетах коэффициент релаксации Сге1х принимался равным 0.5. При вычислении теплового потока со стенок сначала по формулам (2.8) или (2.10) находится Ґ {и на текущей итерации к этому моменту уже известно), после чего значение qw определяется из соотношения, которым вводится величина Ґ.
В этом случае также применяется релаксация с привлечением величины теплового потока с предыдущей итерации. Приведенный алгоритм пригоден и для вычисления массового потока со стенок при решении концентрационного уравнения (2.4). При использовании в качестве замыкающей модели метода LES модели одного уравнения для кинетической энергии неразрешаемого движения, значения ksgs в ближайших к стенкам узлах расчетной сетки находятся из условия (2.6) по известным значениям и т. опыта по применению гибридных подходов к расчету термоконвективных течений была реализована и прошла апробацию простейшая алгебраическая модель, основанная на модели Смагоринского в области LES и соотношениях теории пути смешения Прандтля в области RANS: Как видно, "переключение" в данной модели основано на сравнении вычисляемых в рамках LES и RANS величин турбулентной вязкости. Возможности модели были проанализированы на примере задачи о конвекции Р-Б во вращающемся слое (раздел 3.2). Значительная часть расчетов в рамках настоящей работы (главы 3 и 4) была выполнена на основе дифференциальной гибридной модели RANS/LES, сформулиро
Общая структура течения. Влияние числа Релея. Глобальная циркуляция в полости
Влияние числа Релея на структуру конвекции в полости иллюстрирует рис. 3.25, где представлены мгновенные поля температуры и скорости в одном из меридиональных сечений цилиндра. Видно, что рост числа Релея сопровождается уменьшением толщины температурных пограничных слоев, образующихся у горизонтальных стенок полости. В целом достаточно плавное изменение температуры в вертикальном направлении при наименьшем числе Релея, равном 10 , сменяется сильным градиентным изменением вблизи стенок и близким к постоянству распределением в центральной части полости при Ra=\0s, чему способствуют все более интенсифицирующиеся процессы турбулентного перемешивания (рис. 3.25а, д). В режиме с Ra= 105 наиболее отчетливо проявляется занимающая большую часть полости крупномасштабная конвективная ячейка, состоящая из циркулирующей жидкости (рис. 3.256). Глобальная конвективная структура во всех режимах сосуществует с рециркуляционными зонами в углах полости. Это возвратное движение, образующееся в местах отрыва крупномасштабного течения от стенок полости, способствует оттеснению изотерм от последних, размывая и утолщая тем самым температурные пограничные слои в этих областях. С ростом числа Релея область вовлеченной в циркуляционное движение жидкости со скоростями до 0.5 Vb становится тоньше: ячейка "прижимается" к стенкам полости. Этот процесс сопровождается уменьшением характерных относительных скоростей движения в центральной части полости. Одновременно происходит хаотизация течения: в нем начинают более отчетливо проявляться вторичные вихревые образования (рис. 3.25е). Изменение характера течения с ростом числа Релея наглядно демонстрируется также рис. 3.26, на котором изображены трехмерные картины изоповерхностей температуры (Г =0.25 и 0.75), с нанесенными на них векторами скорости, а также изоповерхности вертикальной компоненты скорости (w = ±0.25) для двух чисел Ra. Как уже было отмечено, переход от умеренного числа Ra = 106 к 108 сопровождается интенсификацией турбулентного перемешивания. Это приводит к заметному искривлению изоповерхностей температуры.
Представленные изоповерхности w, явно указывающие на формирование крупномасштабной конвективной ячейки, в свою очередь, становятся менее "монолитными", распадаясь на отдельные фрагменты. Мгновенные поля и фрагменты полей скорости и температуры в меридиональном сечении полости для режима с числом Ra = 8.7х108, представленные на рис. 3.27, позволяют судить о качестве разрешения динамически важных вихревых структур и высокоградиентных пристенных областей в рамках подхода RANS/LES на данной сетке. Видно, что пограничный слой от стенки до максимума скорости покры вается 7- 8-ю узлами (рис. 3.27а). Метод RANS/LES позволяет прописать на данной сетке и множественные вторичные вихри, формирующиеся как в углах полости, так и в ее центральной части. Кроме того, видно, что изменения температуры порядка 10% от общего температурного перепада происходят на 8 - 9-ти узлах расчетной сетки в пристенных областях у горизонтальных стенок цилиндра (рис. 3.276). Типичные мгновенные поля безразмерной кинетической энергии неразрешаемого движения ksgs и эффективной вязкости veffl v в меридиональном сечении о о цилиндра для двух значений числа Релея 10 и 8.7x10 показаны на рис. 3.28. Безразмерная кинетическая энергия неразрешаемого движения имеет порядок 10"3 с максимумами в местах наиболее интенсивной завихренности потока - в углах области (см. рис. 3.25е, 3.27а). С ростом числа Релея расширяется диапазон неразрешаемых масштабов турбулентности, что в сопоставляемых вариантах сопровождается примерно трехкратным увеличением характерных значений относительной подсеточной вязкости в центральной части полости (рис. 3.286, е). На рис. 3.29 представлены мгновенные поля скорости и ее вертикальной компоненты в среднем горизонтальном сечении полости для трех режимов конвекции. Видно, что во всех случаях сечение приблизительно поровну поделено на области с противоположным знаком вертикальной скорости. С ростом числа Ra граница раздела, соответствующая нулевой скорости вертикального движения, утрачивает четкие очертания, становясь сильно размытой (рис. 3.29д), при этом структура горизонтального движения выглядит все менее упорядоченной (рис. 3.29е). В развитых режимах турбулентной конвекции идентификация крупномасштабной конвективной ячейки по мгновенным картинам поля скорости становится затруднительной (рис. 3.25г, 3.26г). Для преодоления этой проблемы следует исключить из рассмотрения высокочастотные пульсации, фильтрацию которых можно осуществить посредством "пакетного" временного осреднения. Рис. 3.30 иллюстрирует распределения вертикальной компоненты скорости в центральном горизонтальном сечении полости, каждое из которых получено в результате последовательного осреднения ста мгновенных полей (/?а=108). Теперь глобальная структура более отчетливо проявляется, при этом осесимметричная форма полости не препятствует изменению ее азимутального положения. Схематическое изображение пространственной структуры глобальной циркуляции ртути и поля скорости в цилиндрической полости приведено на рис. 3.31.
Конвективная ячейка с линейным размером порядка высоты (диаметра) полости поддерживается движением вдоль стенок когерентных структур (рис. 3.31а). Порождение всплывающих структур становится результатом нагревания жидкости при протекании по горячей области вблизи нижней стенки полости. Обратные процессы наблюдаются у противолежащей холодной стенки. Схема на рис. 3.316 показывает, что крупномасштабное течение имеет погранслойный характер. Приведенные на рис. 3.32а, б профили осредненной во времени максимальной горизонтальной разности температур AT {z) для случаев Ra= 10 и 10 показывают градиентный рост последней в пристенных областях и практически постоянную вели о чину в ядре потока. Однако при Ra 10 профили становятся немонотонными, с плавным уменьшением ATh в направлении центра полости и острыми максимумами в пристенных областях; координаты максимумов смещаются ближе к стенкам с ростом числа Релея. Величина максимальной горизонтальной разности температур в центре полости AThc (z = 0.5) сохраняется в режимах с умеренными числами Релея и слабо уменьшается с ростом числа Релея в соответствие с закономерностью, близкой к Ra 5 при числах Ra 10 (см. рис. 3.33а). Полученные в настоящих расчетах значения AThc находится в хорошем количественном согласии с данными DNS [Verzicco, 1999], где при Ra 106 - AThc 0.6. Оценка характерной скорости глобальной циркуляции Vg может быть осуществлена, например, на основе анализа профиля осредненной во времени максимальной горизонтальной разности вертикальной компоненты скорости Awm. При этом, следуя [Verzicco, 1999] и схеме на рис. 3.316 принимается, что величина Awm в центральном сечении приблизительно равна максимальной вертикальной разности горизонтальной компоненты скорости на оси цилиндра и составляет двукратную величину скорости Vg. Показанные на рис. 3.32в, г профили обнаруживают монотонный рост Awm в направлении центра полости лишь для режима с наименьшим значением числа Релея, а именно при Ra= 105. При дальнейшем развитии турбулентной конвекции на профилях в пристенных областях возникают вторичные максимумы. Аналогично максимальной горизонтальной разности температур (АТ величина Awm в центре полости практически не зависит от числа Релея в режимах с умеренными его значениями и слабо уменьшается пропорционально Ra 005 при Ra 108 (см. рис. З.ЗЗб). Построенные по значениям скорости Vg = 0.5Awmc и высоте полости числа Рейнольдса Reg в зависимости от Ra нанесены на рис. 3.34. Видно, что данные настоящих расчетов достаточно хорошо ложатся на аппроксимационные прямые RegccRa0A5 и согласуются с экспериментальными оценками [Takeshita, 1996; Cioni, 1997; Naert, 1997], произведенными на основе анализа пиковой частоты , на спектре пульсаций температуры (см. раздел 3.1). Иллюстративный график временной эволюции безразмерной максимальной вертикальной компоненты скорости в полости wmax для трех различных Ra приведен
