CC BY
58
УДК 533.6.011:621.165
А.В. Григорьев, А.И. Якунин, Н.Б. Кузнецов, В.Ф. Кондратьев, Н.Н. Кортиков
расчет нестационарного ротор-статорного взаимодеиствия
в турбинной ступени методом гармонического баланса
Пульсации параметров течения, вызванные ротор-статорным взаимодействием в ступенях осевой турбины («clocking» эффект), составляют от 5 до 20 % осредненных по времени значений при типичных осевых зазорах. Столь значительный уровень колебаний скорости потока приводит к тому, что реальные газодинамические характеристики ступени (потери полного давления, КПД, излучаемый тональный шум) отличаются от тех, которые ожидаются в предположении стационарности течения в каждом венце. Изучение пульсаций давления и скорости, связанных с взаимодействием лопаточных венцов, позволяет улучшать аэродинамическое совершенство и надежность будущих поколений авиадвигателей [1].
Пространственное CFD-моделирование лопаточных венцов турбины — самый важный этап в аэродинамическом совершенствовании современных турбин. С развитием вычислительных методов и средств вычислительной техники развивалось моделирование лопаточных венцов — от стационарного решения в одном лопаточном венце до нестационарного решения в многоступенчатых турбомашинах.
Стационарное решение в одном лопаточном венце получается быстро, но оно требует предварительной информации о смежных решетках. Стационарное решение для многоступенчатых турбомашин использует подход «mixing plane», позволяющий связывать между собой смежные лопаточные венцы. Ограничение в данном методе вызвано осреднением параметров потоков в окружном направлении между лопаточными венцами. Несмотря на такие ограничения, расчет многоступенчатых турбомашин с использованием «mixing plane» представляет собой разумный компромисс между точностью и вычислительными ресурсами и, следовательно, остается ценным и мощным инструментом для проектировщика [2].
Поскольку турбины достаточно эффективны, то дальнейшее повышение их эффективности и надежности требует тонкого понимания таких наблюдаемых физических явлений, как распространение кромочного следа, взаимодействия кромочного следа с другими лопатками, ламинарно-турбулентный переход и т. д. Конструктор может использовать разнообразные приемы по улучшению проекта за счет выбора оптимального распределения нагрузки, расположения лопаток в окружном направлении и пр. Однако эффект от каждого способа должен быть оценен в контексте его взаимодействия с нестационарным потоком.
Моделирование нестационарного ротор-ста-торного (transient rotor/stator — TRS) взаимодействия связано с расчетом на сетке, учитывающей изменение во времени положения ротора относительно статора (скользящая сетка — sliding mesh). При этом подходе численное моделирование проводится применительно к полному рабочему колесу или кратному соотношению числа лопаток для статорной и роторной решеток ступени с использованием условий периодичности.
В случае нарушения условия кратности числа лопаток в ступени нестационарные вычисления, в которых возможно применение условия периодичности, выполняются на основе недавно развитого семейства методов расчета в лопаточных венцах, известных как «методы преобразования». В их число входят: «преобразование профиля», «преобразование времени» и «Фурье-преобразование» [3].
В методе «преобразование профиля (PT)», который является частью TRS, коррекция параметров на внутренних границах выполняется автоматически для случая, когда количество лопаток статора не кратно количеству лопаток ротора. В этом приближенном методе используется скользящая сетка, и течение в лопаточных
венцах с различным количеством лопаток моделируется без геометрического масштабирования или изменения геометрии лопатки. Условие периодичности применяется для каждого временного шага, и распределение параметров потока автоматически корректируется поперек ротор-статорных интерфейсов посредством их вытягивания или сжатия для сохранения потоковых величин.
Метод «преобразование времени (ТТ)» может рассматриваться как коррекция метода РТ. Согласно этому методу уравнения движения преобразованы во времени так, чтобы гарантировать периодичность границ. Метод поддерживает точность решения и правильность передачи данных между ротором и статором в случае некратного соотношения лопаток.
В методе «Фурье преобразования (РТ)» история потока на границах расчетной области сохраняется. Для этих целей используют разложение в ряд Фурье, содержащий слагаемые с частотой вращения ротора, а также более высокие гармоники. Этот метод может быть расширен на многоступенчатые турбомашины.
Метод «гармонический баланса (НВ)» был предложен, чтобы учесть динамические нелинейные эффекты и одновременно уменьшить время получения численного решения [4, 5]. Предполагается, что параметры нестационарного газового потока — периодические во времени и пространстве, как это имеет место в тур-бомашинах. Газодинамические параметры ищутся в виде временного ряда Фурье с частотами, кратность которых первоначальной частоте возбуждения выражается целыми числами [6, 7].
Рис. 1. Схема течения в экспериментальной воздушной турбине
В численном алгоритме в качестве искомых переменных выступают коэффициенты ряда Фурье. Ряд Фурье подставляют в нестационарные уравнения газовой динамики, а затем производят их «расщепление». В результате определяют набор упрощенных уравнений, имеющих одну и ту же гармонику. Решение этих уравнений хранится во множестве, соответствующем временным подуровням в пределах одного периода [8, 9].
Цель настоящей работы — апробация метода гармонического баланса (HB) для расчета нестационарных процессов. Апробация основана на сопоставлении результатов расчета с опытными данными для одноступенчатой осевой турбины ЛПИ [10] и данными, полученными по методу TRS с применением скользящей сетки. Все вычисления выполнены с помощью кода STAR CCM + 7.04.006.
Физическое моделирование
Физическое моделирование нестационарных явлений в турбинной ступени проводилось на экспериментальной воздушной турбине ЛПИ ЭТН-4 [10]. При создании стенда особое внимание было уделено обеспечению равномерного по окружности подвода воздуха с малой турбулентностью. Неравномерность потока по окружности перед направляющим аппаратом составляет 1 %, а степень турбулентности — 1,5 %. При создании стенда с целью уменьшения уровня вибраций применялись подшипники качения с увеличенным числом шариков.
Проточная часть турбины выполнена с постоянной высотой канала в меридиональной плоскости. Средний диаметр рабочего колеса — 430 мм при высоте лопаток 70 мм. Схема течения в этой экспериментальной воздушной турбине изображена на рис. 1. Геометрические характеристики решеток приведены в таблице.
Основным источником нестационарности потока в рабочем колесе турбины является неравномерность в кромочных следах за направляющими лопатками. Поэтому большой практический интерес представляет исследование влияния уровня неравномерности в следах на величину переменных аэродинамических сил и на характер нестационарного течения в межлопаточных каналах рабочего колеса.
Рабочие лопатки закреплены на массивном диске. Конструкция диска позволяет устанавли-
Геометрия ступени
Сечение Л, gl, g2, аЬ Р1, Р2,
м м м м м градус градус градус
Корневое 0,36 0,0471 0,0269 0,078 0,035 25°10' 34°16' 36°30'
Среднее 0,43 0,0562 0,0322 0,078 0,037 26°20' 58°25' 34°32'
Периферийное 0,50 0,0654 0,0374 0,078 0,043 28°20' 94°50' 30°16'
Здесь: cí, ^ — хорды сопловой и рабочей лопаток; gí, g2 — шаги в решетках; d — диаметр; а! — угол входа в направляющую решетку; р1; р2 — углы входа и выхода в рабочем колесе.
вать на нем специальные измерительные лопатки для изучения переменных аэродинамических сил. Измеряемая физическая величина преобразуется в электрический сигнал, который передается на вход измерительной и регистрирующей аппаратуры по проводам, уложенным в полом валу ротора, через охлаждаемый водой ртутный токосъемник. На валу ротора установлен диск датчика, вырабатывающего импульсное напряжение с заданной частотой. Датчик построен на основе фотодиода ФД-2.
В опытной турбине обеспечены широкие возможности для измерения шаговой неравномерности и параметров нестационарного потока в четырех сечениях проточной части для 18 измерительных точек (см. рис. 1). С этой целью наружный корпус с закрепленными на нем координатниками мог поворачиваться электродвигателем на 360°вокруг оси турбины. Количе-
ство направляющих лопаток — 24 шт. Количество рабочих лопаток — 48 шт.
Газовый поток в турбине — периодически нестационарный, поэтому необходимо измерение мгновенной скорости. Рис. 2, a показывает форму волны мгновенной скорости в осевом зазоре турбины. Запись формы волны осуществлялась электронным осциллографом, который связан с термоанемометром. В верхней части осциллограммы дается запись таймера сигнала положения лопатки.
Методика обработки осциллограмм заключалась в следующем. Сначала проводилась дискретизация осциллограммы и графическим интегрированием определялась средняя линия, которая соответствует средней скорости. Затем от нее проводился отсчет (в миллиметрах на экране осциллографа) с целью определения величины пульсаций. Измерение переменного во
Рис. 2. Экспериментальные данные для мгновенной скорости в осевом зазоре турбины (а) и положение контрольных точек при численном моделировании (б)
времени давления на поверхности лопатки осуществлялось с помощью датчика давления. Датчик изготовлен в виде гибкой пластины, защемленной по периметру. Величина деформации пластины соответствовала измеряемому переменному давлению.
В процессе эксперимента на выходе рабочего колеса неподвижный зонд термоанемометра делал запись колебаний скорости, которые вызваны следующим: 1) неравномерным потенциальным полем лопаточного аппарата;
2) кромочными следами от рабочих лопаток;
3) возмущениями скорости, возникающими периодически на входе в каждый движущийся межлопаточный канал и проходящими через него, соответствующим образом трансформируясь. Как показано на рис. 2, б, кромочные следы от направляющего аппарата входят в межлопаточные каналы рабочего колеса в виде областей повышенной турбулентной вязкости.
Уровень флуктуаций скорости в точках мониторинга (это, в частности, точки 15, 25, 26, 27 и 28 на рис. 2, б) зависит от предистории, которая определяется особенностями течения газа около вогнутой и выпуклой поверхностей рабочей лопатки, а также геометрией ее выходной кромки. Как можно заметить (рис. 3), левый склон импульса сформирован регулярными
пульсациями, что указывает на стабильный и устойчивый характер течения воздуха около выпуклой поверхности. На рисунке 3 также можно отметить, что уровень колебаний скорости уменьшается с увеличением расстояния от выходной кромки рабочей лопатки.
Значительный уровень представленных на осциллограммах (рис. 3) случайных (нерегулярных) колебаний скорости на правом склоне импульса указывает на высокую интенсивность турбулентности, развивающейся с вогнутой стороны рабочей лопатки.
Этот эффект вызван дестабилизирующим влиянием сил Кориолиса, в результате чего частица жидкости движется в направлении по нормали к вогнутой поверхности лопатки и тем самым увеличивает уровень турбулентных пульсаций около нее. На выпуклой поверхности влияние силы Кориолиса носит стабилизирующий характер.
Вычислительный метод
Нестационарное течение вязкого газа в лопаточных венцах турбомашин описывается уравнениями Навье — Стокса и энергии, записанными в интегральной форме для движущегося контрольного объема в относительной системе координат, вращающейся с угловой скоростью О.
Рис. 3. Экспериментальные данные для мгновенной скорости позади рабочего колеса
Поскольку функция решения W— периодическая во времени, то она может быть представлена рядом Фурье
_ м _ W(г,т) = ^ Wm(г)ехр(/юшт), (1)
m=-M
где ю — фундаментальная частота; М — число гармоник, сохраненных в решении; Wm — коэффициенты ряда Фурье.
Длины и частоты возмущений внутри венцов многоступенчатой турбомашины зависят от количества лопаток в решетках и скорости вращения ротора. В частности, возмущения, генерируемые в результате взаимодействия двух решеток с количеством лопаток Bl и В2, характеризуются окружными волновыми числами (нодальные диаметры), равными N = шфх + + ш2В2, где т1 и ш2 — целые числа. Частоты этих возмущений вычисляются по формуле ю = п2В2О. для неподвижной области, ю = п^О — для вращающейся области.
Дискретные аналоги уравнений для определения коэффициентов ряда (1) выводятся посредством применения метода контрольного объема, разбиения каждого объема на полигек-саэдральные ячейки на основе схемы второго порядка точности. Дискретизация конвективных слагаемых уравнений проводится согласно противопоточной процедуре, а слагаемых, ответственных за перенос диффузией, — центральными разностями. Введение дополнительного слагаемого с локальной производной от псевдовремени позволило использовать маршевый метод и процедуру установления. Решение на
каждом временном подуровне осуществлялось с помощью неявного метода Эйлера, и полученная система линейных алгебраических уравнений решалась с помощью многосеточного AMG-метода.
Условия периодичности в окружном направлении используются с целью уменьшения размеров вычислительной области и сведения к одному или двум межлопаточныи каналам.
Применение «неотражающих» условий на входной и выходной границах расчетной области позволяет избежать неустойчивых численных решений, которые могут возникнуть при наличии случайных возмущений в процессе отражения.
Чтобы получить решение для нулевой моды, соответствующее стационарному решению, применяются традиционные граничные условия, а именно: на входе задаются полные давления, температура и угол натекания потока, на выходе — статическое давление. Решения для высших гармоник строится на основе обратного Фурье-преобразования и результатов численного решения семейства уравнений для коэффициентов ряда (1) на каждом временном уровне.
Результаты расчетов
Метод гармонического баланса использует структурированную сетку типа НОН, которая создана для каждого из двух лопаточных венцов, как показано на рис. 4. Входные и выходные границы для каждого из лопаточных венцов соответствуют плоскостям, где известны экспериментальные данные. Расчетная сетка состоит из
Осевая координата, мм
Рис. 4. Расчетная сетка для метода гармонического баланса и моделирования с применением скользящей сетки (а); меридиональное сечение ступени (б)
1,3 млн ячеек. Размер пристеночной ячейки составляет 0,01 мм.
Условия эксперимента соответствовали следующим величинам: полному избыточному давлению на входе 0,14 МПа; полной температуре 350 K; избыточному статическому давлению на выходе 0,0 МПа; угловой скорости вращения ротора, равной 5000 оборотов в минуту.
Чтобы обеспечить хорошее начальное приближение для ДВ-решателя используется специальная инициализация, в основе которой лежит решение стационарных уравнений Навье — Стокса с интерфейсом в виде «mixing plane», начиная с грубой и заканчивая подробной сеткой. В качестве уравнения состояния принимается уравнение идеального газа; модель турбулентности — модель Спаларта—Аллмареса. Решение осуществлялось с числом Куранта равным пяти в случае одной, трех и пяти мод. Расчет сходился к периодическому нестационарному решению за 2000 итераций.
Рис. 5, а, в показывают картину контуров турбулентной вязкости в среднем сечении лопа-
ток в текущии момент времени, которая получена двумя разными методами (НВ и ТЯ$). Эти рисунки демонстрируют, что кромочный след от направляющей решетки, проходя вращающийся межлопаточный канал, расширяется в продольном и поперечном направлениях.
Этот процесс обусловлен «разломом» кромочного вихря о входную кромку лопатки рабочего колеса и, как следствие, генерацией нестационарных возмущений вверх и вниз по течению. Область потока, расположенная вокруг рабочей лопатки, характеризуется повышенной турбулентной вязкостью, особенно на расстоянии, приблизительно равном половине хорды. Это связано с образованием здесь устойчивой рециркуляционной зоны. Рис. 5, б демонстрирует изменение модуля скорости вдоль границы интерфейса, расположенной в осевом зазоре.
Профиль скорости имеет два минимума: один из них (с большим провалом) ассоциируется с областью следа за направляющей решеткой, а второй — с тормозящим воздействием входной кромки рабочего колеса. На рис. 5, г
б)
V, м/с 210 205
195
190
185
тяз-_
V,
/
4
20
30
40
50
60
70
Y, мм
в)
Рис. 5. Мгновенные картины контуров турбулентной вязкости в срединном сечении лопаток
ступени, полученные с помощью НВ (а) и ТЯ8 (в) методов; распределение скорости по ширине межлопаточного канала внутри осевого зазора (б) и на выходе рабочего колеса (г)
дано распределение скорости вдоль прямой СБ (см. рис. 5, в), которое также содержит два минимума, чье появление связано с кромочными следами от соседних рабочих лопаток. Потенциальная неравномерностть потока в рабочей решетке и деформированный кромочный след от направляющего аппарата определяют немонотонное изменение профиля скорости между этими двумя провалами. Рис. 6, а, б, в, г показ-вают распределение осредненного во времени (нулевая мода) давления и поправки на нестационарность для давления (первая мода) по обводу лопаток статора и ротора одноступенчатой турбины ЛПИ, подсчитанные по методу гармонического баланса.
Результаты расчетов показывают хорошое согласование с опытными данными для осред-ненного во времени давления по обводу статор-ной лопатки (рис. 6, а). Результаты расчетов, представленные на рис. 6, б, в, иллюстрируют малое влияние нестационарности на величину нулевой моды давления на поверхности ротор-
ной лопатки, а также первой моды на поверхности статорной. Важно отметить, что пульсации давления (изменение величины первой моды) на поверхности роторной лопатки могут достигать (вблизи входной кромки лопатки рабочего колеса — рис. 6, г) восьми процентов от осредненного значения давления.
Сравнение результатов вычисления с помощью TRS и HB подходов дано на рис. 7—9. Видно, что различие между двумя методами для мгновенной скорости незначительно (меньше пяти процентов), если число мод в методе гармонического баланса выбрано равным пяти.
Сравнение результатов расчета с экспериментальными данными осциллограмм (для точек мониторинга 76 и 28) представлено на рис. 8 и 9. Оба метода (TSR и HB) удовлетворительно прогнозируют форму волны в поведении модуля скорости с течением времени, отмечая отличие в поведении кривых от чисто гармонической функции. Но оба метода не могут отследить
а)
Давление на поверхности, Па
-0,04 -0,03 -0,02 -0,01 0 0,01 0,02 Осевая
координата, м
б)
Давление на поверхности, Па
^Г-....... N
\
—---
0,04 0,05 0,06 0,07 0,08 Осевая
координата, м
в)
Давление на поверхности, Па
Real ч . А
, %
A-J : Ч*
/ i
Imag, /
■j
-0,04 -0,03 -0,02 -0,01 0 0,01 0,02 Осевая
координата, м
г)
Давление на поверхности, Па
GBO 100 iflO 0.0 -2 ае-
-4»»
Real У
- ......У
wr> Л
ь У 7
■ \ Vх
\ Imag
V /
0,04 0,05 0,0
0,07 0,08 Осевая
координата, м
Рис. 6. Нулевая мода в Фурье-представлении давления на поверхностях лопатки статора (а) и ротора (б); реальные и мнимые части первой моды для статора (в) и ротора (г)
а)
V, м/с
110 105 100
95 90
TRS /
\ /у х/ \
-Л /
V 7 НВ \7
V 1 V
б) V, м/с
115
110
105
100
95
90
НВ TRS
/
/
/
0.0473 0.0474 0.0475 0.0476 0.0477 0.0478 Время, с
¿.а;- .■'.¿■Г4 .¿.'М-':. ¿¿-.Г'" 0.04Г-- ¿.¿;-з Время, с
а)
V, м/с
190
185
180
Рис. 7. Сравнение результатов вычислений для TRS и HB подходов для двух точек рис 5, а:
15 (а) и 27 (б)
б)
/ \ НВ
/ \ // TRS \ 1
/
0,0473 0,0474 0,0475 0,0476 0,0477 Время, с
Рис. 8. Зависимость скорости от времени на расстоянии 8,0 мм от выходной кромки статора для точки 76 на рис. 1 (а — моделирование; б — экспериментальные данные)
а)
V, м/с 110 105 100 95
б)
90
0,0473 0,0474 0,0475 0,0476 0,0477 Время, с
Рис. 9. Зависимость скорости от времени на расстоянии 10,0 мм от выходной кромки рабочей лопатки для точки 28 на рис. 1 (а — моделирование; б — экспериментальные данные)
эффекты нестационарности во всем диапазоне частот. В частности, на расчетных кривых для скорости потока отсутствуют высокочастотные участки, характерные для малых масштабов турбулентности.
Метод гармонического баланса (НВ), встроенный в пакет STAR-CCM +, применен для расчета нестационарных эффектов при течении вязкого газа в одноступенчатой турбине ЛПИ. Показана приемлемая точность НВ-метола для расчета регулярных колебаний потока в осевом зазоре и на выходе рабочего колеса. Такой вывод сделан на основе сравнения результатов расчета
с опытными данными и результатами численного моделирования, в котором применялась полноразмерная модель ступени и скользящих сеток: погрешность не превышает пяти процентов. Отмечается, что пульсации давления на поверхности роторной лопатки в окрестности входной кромки могут достигать восьми процентов относительно осредненного во времени значения давления. Дается объяснение нелинейных эффектов в поведении профиля скорости в межлопаточных каналах ступени, связанных с ротор-статорным взаимодействием и развитием кромочных следов за направляющей решеткой и рабочим колесом.
СПИСОК ЛИТЕРАТУРЫ
1. Breugelmans, F. Resent research in the VKI tur-bomachinery and propulsion department [Текст] / F. Breugelmans // Сб. тез. Второй международной научно-техн. конф. «Авиадвигатели XXI века».— Москва.— 2005. Т. 1.— С. 16-25.
2. Biesinger, T. Unsteady CFD methods in commercial solver for turbomachinery applications [Текст] / T. Biesinger, Ch. Cornelius, C. Rube [et all.] // Proc. of ASME Turbo Expo.— Power for Land, Sea and Air. Glasgow, 2010.— UK. 12 p.
3. Connell, S. A comparison of advanced numerical techniques to model transient flow in turbomachinery blade rows [Текст] / S. Connell, M. Braaten, L. Zori [et all.] // Proceedings of ASME Turbo Expo— 2011. Vancouver, British Columbia, Canada.— 2011. 10 p.
4. He, L. Efficient approach for analysis of unsteady viscous flows in turbomachines [Текст] / L. He, W. Ning // AIAA J.— 1998. Vol. 36, № 11.— P. 005-2012.
5. He, L. Analysis of rotor— rotor and stator -stator interferences in multi— stage turbomachines [Текст] / L. He, T. Chen, R.G. Wells [et all.] // Proc. of GT2006 ASME Turbo Expo— 2006.— Power for Land, Sea and Air. Barcelona, Spain. 2006.— 11 p.
6. Hall, K.C. Computation of unsteady nonlinear flows in cascades using a harmonic balance technique [Текст] / K.C. Hall, J.P. Thomas, W.S. Clark // AIAA J.— 2002.Vol. 40, № 5.— P. 879-886.
7. Vilmin, S. Unsteady flow modeling across the rotor/ stator interface using the nonlinear harmonic method [Текст] / S. Vilmin, E. Lorrain, Ch. Hirsch, M. Swoboda // Proc. of ASME Turbo— Expo 2009.— Power for Land, Sea and Air. Orlando, Florida, 2009. USA.— 9 p.
8. Weiss, J.M. Simulation of unsteady turbomachinery flows using of implicitly coupled nonlinear harmonic balance method [Текст] / J.M. Weiss, V. Subramanian, K.C. Hall // Proc. of ASME. Turbo Expo— 2011.— Vancouver, Canada, 2011.— 8 p.
9. Custer, C.H. Unsteady simulation of a 1.5 stage turbine using an implicitly coupled nonlinear harmonic balance method [Текст] / C.H. Custer, J.M. Weiss, V. Subramanian [et all.] // Proc. of ASME Turbo Expo— 2012.— Copenhagen. Denmark, 2012.— 15 p.
10. Кондратьев, В.Ф. Исследование физических явлений и нестационарных взаимодействий в турбинной ступени [Текст]: Автореф. ... канд. техн. наук / В.Ф. Кондратьев.— Л., 1977.— 18 с.
