Теплофизика высоких температур, 2023, T. 61, № 1, стр. 75-82

ВВЕДЕНИЕ

В настоящее время подавляющая доля электроэнергии вырабатывается на тепловых и атомных электростанциях с помощью турбогенераторов. При работе турбогенератора в обмотках ротора, статора и активной стали выделяется тепло, которое необходимо от них отводить. Превышение температуры активных частей над классом нагревостойкости изоляции приводит к ее ускоренному тепловому старению, что может послужить причиной электрического пробоя и привести к возникновению аварийных ситуаций на электростанциях. Поэтому вопросы разработки, проектирования и оптимизации системы охлаждения турбогенератора очень важны. В инженерной практике применяются различные подходы для расчета системы охлаждения, в том числе используются численные решения уравнений Рейнольдса (RANS).

Система охлаждения ротора турбогенератора состоит из множества прямоугольных вращающихся каналов, одним из типичных ее элементов является канал с поворотом на 90°.

Опыт исследований влияния вращения и поворота на течение в канале достаточно обширен. В частности, течения в прямых вращающихся каналах, ориентированных параллельно (аксиальных) и перпендикулярно (радиальных) оси вращения, хорошо изучены численно [1‒10] и экспериментально [11‒13]. В этих каналах исследованы структура потока, влияние вращения на эту структуру и турбулентные характеристики, а также корреляции числа Нуссельта на стенках в зависимости от безразмерных критериев течения. Существует множество работ по исследованию течения в круглых каналах с поворотом на 90° [14‒17], а также во вращающихся каналах с поворотом на 180° [18‒22], где предприняты попытки оптимизации формы таких каналов для снижения их сопротивления. Кроме того, показана возможность создания анизотропных RANS-моделей турбулентности с улучшенными возможностями предсказания влияния вращения.

Примером исследования течения в каналах квадратного сечения с резким поворотом на 90°, где отрыв потока вызван скорее формой канала, нежели неблагоприятным градиентом давления, является работа [23], в которой приведено сравнение данных эксперимента, прямого численного моделирования и RANS. Различие в профилях скорости и турбулентных характеристиках за поворотом составляет до 40%, что говорит о посредственной точности моделирования течений в каналах такого типа с помощью RANS-подхода. Касательно течений во вращающихся каналах прямоугольного сечения с поворотом на 90° авторам не удалось найти каких-либо экспериментальных или численных исследований. Данная работа призвана хотя бы отчасти восполнить этот пробел.

Таким образом, целью работы является оценка точности расчета течения во вращающемся прямоугольном канале с поворотом на 90° методом RANS. К сожалению, натурный эксперимент в роторе вращающегося турбогенератора позволяет получить только интегральные характеристики потока в каналах системы охлаждения, что недостаточно для оценки точности RANS-расчетов. Поэтому в данной работе проведен расчет существенно более точным по сравнению с RANS вихреразрешающим методом WMLES, результаты которого использованы в качестве “эталонных”.

ПОСТАНОВКА ЗАДАЧИ И МЕТОД РЕШЕНИЯ

Размеры канала и параметры рассматриваемого течения выбраны таким образом, чтобы соответствовать средним значениям в каналах ротора типичного турбогенератора. В частности, выбранные размеры сечения канала соответствуют сечению вентиляционных каналов в роторе турбогенератора с воздушным охлаждением средней мощности, а угловая скорость вращения канала равна угловой скорости вращения такого ротора.

Рассматривается течение воздуха в канале с поворотом на 90°, вращающемся вокруг оси оХ (рис. 1), расположенной в плоскости симметрии канала. Расстояние от оси вращения до центра аксиального (ориентированного вдоль оси оХ) канала составляет R. Прямоугольное сечение канала имеет размеры H × 2.78H, где H = 0.125R. В роторе реального турбогенератора длина канала перед поворотом достаточна для установления развитого режима течения. Длина аксиальной части канала составляет 3H, и течение в окрестности входной границы полагается установившимся. Выходная граница канала расположена на расстоянии 30H от поворота и имеет неплоскую форму: все точки ее поверхности находятся на одном радиусе от оси вращения.

Рис. 1.

Вращающийся канал прямоугольного сечения с поворотом на 90°.

Газ полагается совершенным с постоянными вязкостью и теплопроводностью, не зависящими от температуры. Число Рейнольдса, посчитанное по среднемассовой скорости V_bulk, высоте канала H и кинематической вязкости газа ν, равно Re = HV_bulk/ν = 14 500 (что примерно соответствует числу Рейнольдса, построенному по динамической скорости, Re_τ = 415), число Маха на входе в канал построено по среднемассовой скорости и скорости звука на входной границе a: M = V_bulk/a = 0.12 (течение слабосжимаемо).

Для оценки влияния эффектов вращения в работе рассматривается как течение в неподвижном канале, так и течение в канале, вращающемся с угловой скоростью ω = 100π рад/c, чему соответствует число Россби Ro = ωH/V_bulk = 0.035. Для течений с таким числом Россби характерны явления стабилизации и дестабилизации течения на сторонах пониженного и повышенного давления соответственно, а также возникновение вихрей Тейлора‒Гетлера [24].

Температура воздуха на входе в канал составляет T_b = 340 К, температура стенок постоянна и равна T_w = 400 К. Число Грасгофа, построенное по центробежной силе, равно Gr = β∆Tω²RH³/ν² = = 4.5 × 10⁸. Число Ричардсона Ri = Gr/Re² ≈ 2.14 близко к единице, что свидетельствует о смешанном режиме конвекции [25]. Для учета влияния свободной конвекции рассматривается именно сжимаемая жидкость.

Так как целью данной работы является оценка точности RANS-моделей турбулентности в сравнении с вихреразрешающим расчетом, задача решалась в двух постановках. В первой из них решаются стационарные уравнения Рейнольдса, замкнутые с использованием различных моделей турбулентности. Во второй постановке использовался вихреразрешающий подход – LES с пристеночным RANS-моделированием (WMLES).

Для моделирования турбулентности в RANS-расчетах использовались модели турбулентности, которые способны учитывать влияние вращения на поток: SA RC [26], k–ε RNG RC EWT [27], k–ω SST RC [28] и BSL-EARSM [29]. WMLES [30] предполагает решение отфильтрованных уравнений Навье‒Стокса, замкнутых с использованием подсеточной модели Смагоринского в большей части расчетной области (вне пограничного слоя и во внешней его области), в оставшейся части пограничного слоя для замыкания уравнений движения применяется алгебраическая модель Прандтля.

Уравнения движения рассматриваются во вращающейся неинерциальной системе координат, поэтому в них в явном виде присутствуют центробежная и кориолисова силы. Поскольку сила тяжести много меньше центробежной силы, она не рассматривается.

На входной границе задаются профиль скорости и турбулентные характеристики развитого течения, посчитанные по модели BSL-EARSM в несжимаемой постановке, и постоянная температура T_b = 340 К. В WMLES турбулентные пульсации на входе создаются с помощью генератора синтетической турбулентности STG [31]. На выходной границе задается постоянное давление. Строго говоря, в общем случае давление на выходной границе канала не является постоянным из-за действия центробежной и кориолисовой сил. Благодаря цилиндрической форме выхода влияние первой из них отсутствует, влиянием силы Кориолиса пренебрегается. Дополнительные расчеты показали, что это допущение влияет на течение лишь в малой окрестности выхода (<1H) и не затрагивает область течения, представляющую интерес. На твердых стенках задаются условие прилипания и постоянная температура T_w.

Симметричные относительно плоскости симметрии канала многоблочные структурированные расчетные сетки с первым пристенным шагом h₁ = = 10^–3H обеспечивают безразмерное расстояние от стенки до центра первой ячейки y⁺ < 1 по всей поверхности стенок за исключением малой окрестности острой кромки поворота. Шаги расчетной сетки в центральной части канала для WMLES удовлетворяют требованиям для подобных расчетов: ∆x = δ/10, ∆z = δ/20, где Δz ≤ Δy ≤ Δx ‒ максимальные размеры ячейки в направлениях x, y, z соответственно, а δ = H/2 – толщина пограничного слоя в развитом течении. Общий размер сеток для RANS- и WMLES-расчетов составил 1.1 и 5.2 млн ячеек соответственно.

Задача решалась в прикладном пакете ANSYS Fluent 19R1. В RANS-расчетах для решения системы уравнений использовался сопряженный решатель coupled [32]. Пространственная дискретизация уравнений производилась по противопоточной схеме второго порядка SOU [33]. Для интерполяции давления на грань использовалась схема PRESTO! [34]. Интерполяция плотности и величин, замыкающих уравнения Рейнольдса, производилась по схеме SOU. Для вычисления градиентов применялась схема Грина‒Гаусса GGCB [35], записанная относительно центров ячеек. В случае наличия вращения в процессе решения для улучшения сходимости постепенно уменьшались число Куранта и коэффициенты релаксации для давления и количества движения.

В WMLES-расчетах для пространственной дискретизации уравнений сохранения момента импульса и энергии использовались центрально-разностные схемы (CD), для интерполяции давления ‒ схема PRESTO!, а для градиентов – схема GGCB.

В качестве начального приближения использовалось решение, полученное по модели BSL-EARSM с наложенными с помощью внутренних средств Fluent турбулентными пульсациями. Шаг по времени выбирался исходя из выполнения условия ограничения числа Куранта единицей в большей части расчетной области (в зонах сгущения сетки число Куранта незначительно превышает единицу) и составил Δt = 0.016T_conv, где T_conv = H/V_bulk.

Осреднение проводилось после установления средних значений среднемассовой скорости и трения на стенках на протяжении T_ave = 185T_conv для течения без вращения и с вращением.

Расчеты проводились на вычислительном кластере Санкт-Петербургского политехнического университета “Торнадо” с использованием двух узлов, состоящих из 56 процессоров Intel Xeon с тактовой частотой 2.6 ГГц каждый, и занимали порядка 100 ч на один вихреразрешающий расчет.

СТРУКТУРА ТЕЧЕНИЯ И ВЛИЯНИЕ ВРАЩЕНИЯ

Структура течения, полученная по результатам WMLES-расчетов для неподвижного и вращающегося каналов, представлена на рис. 2. В неподвижном аксиальном канале (рис. 2а, цифра 1) в развитом течении присутствуют четыре основных и угловых вторичных вихря, практически симметричных относительно плоскостей симметрии аксиального канала. Часть потока отрывается от нижней стенки в окрестности внешнего угла поворота, при этом возникает вихрь 2. При повороте поток отрывается от внутреннего угла поворота, и образуется отрывная зона, состоящая из двух симметричных вихрей 3, размер которой можно оценить по пристеночным линиям тока на передней стенке канала (причиной неидеальности симметрии является конечное время осреднения решения). После поворота структура вторичных токов меняется, образуется два основных вихря 4, 5, аналогичных вихрям Дина, возникающим в круглых каналах с поворотом.

Рис. 2.

Общая структура течения без вращения (а) и с вращением (б): I – входная граница, II – сечение Z = 0, III – сечение Y = R + 5.5H, IV – сечение Y = R + 25.5H.

Основное влияние вращения выражается в существенном нарушении симметрии (рис. 2б). Отрыв 3 становится более протяженным у стороны пониженного давления и менее протяженным у противоположной стороны, однако площадь отрыва, определяемая по пристеночным линиям тока на передней стенке, не изменяется. После поворота вращающееся течение устанавливается дольше, образуется два крупных вторичных вихря 5, вытянутых вдоль длинной стенки канала. В углах, смежных со стороной повышенного давления, вторичных токов не образуется, размер угловых вторичных токов на стороне пониженного давления 6 меньше, чем в неподвижном канале.

Сравнение осредненного числа Нуссельта, вычисленного по формуле Nu = 〈q_w〉H/(λ(T_w – T_b)), на стенках для случаев без вращения и с вращением свидетельствует о том, что при повороте интенсивность теплоотдачи изменяется при наличии вращения, возникает асимметрия (рис. 3). На стороне повышенного давления число Нуссельта примерно такое же, как и в неподвижном канале, но на противоположной стороне интенсивность теплоотдачи с удалением от поворота падает из-за затухания турбулентных пульсаций скорости и температуры.

Рис. 3.

Распределение числа Нуссельта на стенках неподвижного (а) и вращающегося (б) каналов.

Несмотря на увеличение максимальной длины отрывной зоны, ни ее площадь, ни занимаемый вихрями объем не меняются значительно. Тем не менее перепад давления уменьшается как на участке отрыва и перестроения потока, так и на всем участке поворота (табл. 1, сечения 2‒3 и 1‒3 соответственно). Это связано с тем, что при наличии вращения центробежная сила способствует повороту течения в радиальный канал. Здесь 1 – сечение x = 3H, 2 – сечение y = R + 0.5H, 3 – сечение y = R + 5.5H.

При наличии вращения средняя интенсивность теплоотдачи на стенках уменьшается на 7%, перепад температуры на повороте увеличивается на 15%. Согласно многочисленным исследованиям течений во вращающихся каналах, вращение приводит к значительной интенсификации теплообмена. В данном течении этого не происходит по причине существенных дополнительных вторичных токов на стороне повышенного давления (вихрей Тейлора–Гетлера). Однако наличие вращения приводит к уменьшению пульсаций скорости и температуры и числа Нуссельта соответственно на стороне пониженного давления по сравнению с неподвижным каналом.

ОЦЕНКА ТОЧНОСТИ МОДЕЛЕЙ ТУРБУЛЕНТНОСТИ

Течение в канале с поворотом без вращения. При повороте все рассмотренные полуэмпирические модели турбулентности (SA RC, k–ε RNG RC EWT, k–ω SST RC, BSL-EARSM) существенно завышают протяженность отрыва на передней стенке канала (рис. 4 и табл. 2) по сравнению с результатом WMLES. Это приводит к значительному увеличению перепада давления на участке отрыва и перестроения потока. Однако чувствительность перепада давления на всем участке поворота к длине отрывной зоны невелика (табл. 2).

Рис. 4.

Структура течения без вращения, предсказываемая моделями: (а) ‒ SA RC, (б) ‒ k–ε RNG RC EWT, (в) ‒ k–ω SST RC, (г) ‒ BSL-EARSM.

В радиальной части канала только нелинейная модель турбулентности BSL-EARSM правильно предсказывает наличие угловых вторичных токов, хотя крупные вихревые структуры смещает к передней стенке. Что касается остальных моделей, они правильно предсказывают размер и положение крупных вторичных вихрей.

По сравнению с вихреразрешающим расчетом RANS-модели существенно занижают интенсивность теплоотдачи в области отрыва на передней стенке (рис. 5). Среднее число Нуссельта при повороте (на стенках между сечениями x = 3H и y = R + 5.5H) ниже на 11‒20%, поэтому разница температур меньше (табл. 2).

Рис. 5.

Распределение числа Нуссельта на повороте для моделей: (а) ‒ SA RC, (б) ‒ k–ε RNG RC EWT, (в) ‒ k–ω SST RC, (г) ‒ BSL-EARSM.

Течение в канале с поворотом с вращением. Как уже было показано, при наличии вращения вихреразрешающий расчет предсказывает увеличение максимальной длины отрывной зоны без изменения занимаемой ею площади (рис. 2). Из RANS-моделей турбулентности это воспроизводит только модель турбулентности BSL-EARSM (рис. 6г), линейные модели турбулентности, напротив, предсказывают как уменьшение длины и занимаемой площади отрывной зоны, так и уменьшение перепада давления при повороте по сравнению с неподвижным каналом (табл. 3). Кроме того, линейные модели турбулентности не воспроизводят правильной структуры вторичных токов в радиальной части канала (рис. 6).

Рис. 6.

Структура течения с вращением, предсказываемая моделями: (а) ‒ SA RC, (б) ‒ k–ε RNG RC EWT, (в) ‒ k–ω SST RC, (г) ‒ BSL-EARSM.

Столь различная реакция на вращение приводит к улучшению согласования моделей турбулентности и вихреразрешающего расчета. Распределение чисел Нуссельта при повороте (на стенках между сечениями x = 3H и y = R + 5.5H) с наличием вращения приближается к результатам вихреразрешающего моделирования (рис. 7), для всех моделей турбулентности точность предсказания среднего числа Нуссельта при повороте вырастает на 1‒4%, как и точность предсказания разности температур (табл. 3). Тем не менее нельзя исключать, что такое улучшение согласования является результатом взаимной компенсации ошибок в RANS. С изменением скорости вращения ошибка RANS-моделей может возрасти, а согласование с вихреразрешающим расчетом – уменьшиться.

Рис. 7.

Распределение числа Нуссельта на повороте для моделей: (а) ‒ SA RC, (б) ‒ k–ε RNG RC EWT, (в) ‒ k–ω SST RC, (г) ‒ BSL-EARSM.

ЗАКЛЮЧЕНИЕ

В работе выполнено численное исследование турбулентного течения во вращающемся канале прямоугольного сечения с поворотом на 90°. Проведен численный эксперимент методом WMLES, исследовано влияние вращения на течение. Изучена точность RANS-подхода с использованием различных полуэмпирических моделей турбулентности путем сравнения с результатами численного эксперимента.

Показано, что при наличии вращения вихреразрешающий расчет предсказывает смещение отрывной зоны к стороне пониженного давления и уменьшение потерь давления при повороте. Средняя интенсивность теплоотдачи при повороте уменьшается на 7%, разность температур увеличивается на 15% по сравнению с неподвижным каналом.

Установлено, что при отсутствии вращения независимо от используемой модели турбулентности длина зоны отрыва, предсказанная методом RANS, существенно завышена. Коэффициент теплоотдачи при этом занижается на 11‒21% и разность температур при повороте завышается на 4‒13% по сравнению с результатами WMLES.

Влияние вращения на результаты RANS диаметрально противоположно влиянию, наблюдаемому в WMLES. Полуэмпирические модели турбулентности предсказывают уменьшение отрывной зоны при наличии вращения. Это приводит к улучшению согласования интегральных характеристик при повороте с результатами вихреразрешающего расчета. Возможно, такое улучшение является результатом компенсации ошибок в RANS.

Наименьшую точность при расчете неподвижных и вращающихся течений в прямоугольном канале с поворотом на 90° дает SA RC-модель турбулентности с одним уравнением для замыкания уравнений Рейнольдса. Модели с большим количеством уравнений имеют приблизительно одинаковую предсказательную способность, их точность ограничена скорее самим RANS-подходом, нежели особенностями построения конкретной модели.

Расчеты выполнены с использованием вычислительных ресурсов суперкомпьютерного центра СПбПУ “Политехнический” при финансовой поддержке Минобрнауки России в рамках реализации программы Научного центра мирового уровня по направлению “Передовые цифровые технологии” СПбПУ (соглашение № 075-15-2020-934 от 17.11.2020 г.).