УДК 621.51.226.2.53
А.А. Хориков, В.В. Шкуров, В.А. Шорстов
ФГУП «Центральный институт авиационного моторостроения имени П.И. Баранова»
К ВОПРОСУ О ПРОГНОЗИРОВАНИИ ФЛАТТЕРА ЛОПАТОК БЛИСКОВ КОМПРЕССОРОВ
Предложен метод расчета устойчивости к изгибно-крутильному флаттеру лопаток блисковых конструкций рабочих колес, основанный на вычислении работы нестационарных аэродинамических сил в вязком потоке за период колебаний. Обоснование метода выполнено на основе математической модели флаттера, которая учитывает механическую и аэродинамическую связанность лопаток в колесе. Для расчета нестационарных аэродинамических сил использован разработанный в ЦИАМ программный комплекс Cobra. Валида-ция предложенного метода проведена по экспериментальным данным, полученным при испытаниях вентилятора на стенде ЦИАМ. Учет влияния деформации диска, сдвига фаз между изгибной и крутильной составляющими колебаний и конструкционного демпфирования позволил получить хорошее совпадение с экспериментом.
Ключевые слова: лопатки, блиски, флаттер, частоты и формы колебаний, нестационарные силы, сдвиги фаз, энергетический метод, демпфирование, вязкое течение.
Введение
Большинство известных двигателестроитель-ных фирм при создании новых двигателей для обеспечения отсутствия в эксплуатации флаттера лопаток используют программные комплексы, основанные на статистическом обобщении накопленного экспериментального материала. Однако, объем экспериментальных данных по случаям возникновения флаттера блисков пока весьма незначителен, что не позволяет создать программный комплекс, аналогичный имеющимся для лопаток с замковым соединением. Поэтому в данной работе предлагается метод прогнозирования устойчивости к флаттеру лопаток блисковых конструкций рабочих колес, основанный на вычислении работы нестационарных аэродинамических сил за период колебаний при наличии механической и аэроупругой связанности лопаток в колесе с учетом сдвига фаз между лопатками, а также сдвига фаз между изгибной и крутильной составляющими колебаний.
1. Математическая модель и алгоритм расчета
Кинематика и динамика решетки профилей схематически ранее были представлены в работе [1], где каждый из профилей может совершать как поступательные - л, так и угловые - Э перемещения, обусловленные механическими и аэродинамическими связями. Работа нестационарных аэродинамических сил, при таких изгибно-кру-тильных колебаниях определяется по формуле
где Б и М — векторы нестационарной аэродинамической силы и момента;
Л и - — векторы скоростей перемещения и вращения профиля.
Интегрирование производится за период колебаний профиля. При этом учитываются дополнительные силы и моменты, возникающие из-за несинфазности перемещения и вращения профиля.
Результаты расчетов флаттера, в основном, представляются в виде зависимостей коэффициента аэродинамического демпфирования, который рассчитывается по формуле
^аэр
А
4 • Ек
(2)
где А — работа нестационарного аэродинамического потока;
Ек — кинетическая энергия профиля. Кинетическая энергия для профиля, совершающего изгибные колебания, рассчитывается по формуле
2
Eкизг = mh /2 ,
(3)
где т и Л — масса и скорость профиля соответственно.
Для профиля, совершающего крутильные колебания:
A — J (Fh + MJJ)dt,
(1)
E — T J2 / 2
^ккруч ' ^ ■
(4)
© А.А. Хориков, В.В. Шкуров, В.А. Шорстов, 2012
0
где I и - момент инерции и угловая скорость профиля соответственно.
Для профиля, совершающего изгибно-крутиль-ные колебания:
f = F + F
к к изг к круч'
(5)
Составной частью предлагаемого метода является расчет частот и форм колебаний системы лопатки — диск с бегущими волнами деформации. При этом необходимо учитывать, что при реализации бегущих волн деформации на диске реализуются формы с узловыми диаметрами, которые согласно [1] будут увеличивать крутильную составляющую колебаний лопатки. Авторы отчета для расчета частот и форм колебаний использовали программный комплекс ANSYS Release 13.0 с надстройкой Workbench 2.0 Framework. Результаты использования этого комплекса для расчета частот и форм колебаний демонстрируются на рис. 1. На рис. 2 представлено приращение крутильной составляющей колебаний для профиля на 90% высоты лопатки, влияние которого затем подробно исследуется в следующем разделе. В данном случае приращение составило почти 30%, хотя частота колебаний при этом изменилась всего на 2%.
Рис. 1. Результаты расчета формы колебаний лопаточного венца с диском
Для расчета двумерного вязкого течения в межлопаточных каналах был использован программный комплекс Cobra [2], который хорошо зарекомендовал себя при решении задач внешней и внутренней аэродинамики, а также при решении отдельных задач аэроупругости лопаточных машин. Очевидно, что исходными данными для такого расчета должны быть параметры стационарного течения и геометрии решетки: расход воздуха через рабочее колесо, полное давление, полная температура, густота решетки, количество лопаток в колесе.
В данной постановке не будет учитываться неоднородность лопаток в колесе, которая обычно является стабилизирующим фактором при рассмотрении устойчивости лопаток к флаттеру.
2. Результаты расчетов
2.1. Расчет частот и форм колебаний и его сравнение с экспериментальными данными
Была построена модель сектора диска с лопаткой. К модели в дальнейшем применялись условия поворотной симметрии. Конечно-элементная сетка строилась автоматически инструментами надстройки «ANSYS Workbench», при этом для лопатки использовался алгоритм «Patch Independent» метода «Tetrahedrons», для элемента диска — алгоритм «Patch conforming» того же метода с последующим определением условий контакта по сопрягаемым поверхностям.
Закрепление осуществлялось по цапфе диска в окружном и осевом направлениях. Свойства материала лопаток соответствовали свойствам сплава «Ti64» при температуре 20 °С .
На рис. 4 представлено сравнение результатов расчета с экспериментальными данными, из которого следует, что имеется очень хорошее совпадение.
400 200 0
0 10 20 30 40 50 60 70 80 90 100
Рис. 2. Приращение крутильной составляющей колебаний при учете влияния диска
Рис. 3. Расчетные и экспериментальные значения частот колебаний рабочих лопаток 1-й ступени -- расчет ЦИАМ, ANSYS Release 13.0,
• - экспериментальные значения, стенд ЦИАМ
2.2. Расчет стационарного течения
Начальными данными для исследования обтекания колеблющихся лопаток являются параметры стационарного течения в исследуемой решетке профилей.
Исходные данные для р
Исходными данными для исследования течения в решетке являлись определенные в эксперименте интегральные параметры потока. Их величины и принятая эквивалентная площадь проходного сечения представлены в таблице 1.
Таблица 1
ета параметров течения
Параметр Обозн., ед. измер. Величина
Полное давление на выходе Р*, Па 127530
Полная температура на входе Т*, К 323
Расход воздуха через рабочее колесо Fm, кг/с 55,5
Площадь расчетного канала S, м2 0,346523
Средняя приведенная осевая скорость X 0,376
Число лопаток венца 24
Осевая приведенная скорость потока на входе в решетку 1, необходимая для исследования течения в решетке, определялась по уравнению расхода, площади канала и полным параметрам.
Решетка была составлена из 24-х профилей, являющихся сечением пера лопатки цилиндром радиуса Я = 0,3496 м (~90% высоты лопатки). Полные параметры на входе в решетку
Исходные данные для расчета стаци
берутся из исходных данных для трехмерного течения.
Расчеты стационарного течения производились для двух значений расходов воздуха через венец: основного О1 = 55,5 кг/с (без флаттера) и пониженного О2 = 52,6 кг/с (с флаттером). Значения параметров потока на входе в решетку представлены в таблице 2.
Таблица 2
арного течения в решетке профилей
Параметр Обозн., ед. измер. Величина для 01 Величина для 02
Статическое давление на входе Г, Па 117331 118463
Осевая составляющая абсолютной скорости потока на входе С1а, м/с 123,5855 116,3306
Окружная составляющая абсолютной скорости потока на входе С1и, м/с 326,561 326,561
Плотность р,кг/м3 1,29596 1,30488
Температура Т, К 315,4 316,3
Осевая приведенная скорость 0,376 0,354
2.3. Расчет нестационарного течения
Расчет нестационарного течения в решетке профилей проводился для всех 24-х межлопаточных каналов.
Использовалась подвижная (отслеживающая движение профилей) расчетная сетка, содержащая около 640000 ячеек. Условие обобщенной периодичности не использовалось.
Движение сетки осуществляется в выделенных областях, каждая из которых представляет собой топологический прямоугольник, содержащий исследуемый профиль. Границы четырехугольной области остаются неподвижными (рис. 4).
Рис. 4. Расчетная область
В исследовании было принято, что лопатка колеблется только на одной (первой) собственной частоте. Движение выбранного для расчета профиля по первой форме колебаний соответствует его поступательному перемещению ^ и повороту на некоторый угол д.
В соответствии с применяемым для исследования энергетическим методом анализа флаттера, в задаче моделировались гармонические колебания профилей.
Для начальных расчетов амплитуды были подобраны так, что соотношение между изгибом и кручением профиля соответствовало соотношению, реализующемуся при движении профиля по первой форме колебаний на 90% высоты изолированной лопатки, т.е. без учета влияния диска. Таким образом, было принято ^ = 0,001 м, а д = 0,0025 рад.
Параметрами задачи являлись:
— сдвиг фаз между перемещениями соседних профилей, соответствующими изгибным колебаниям лопатки -
— сдвиг фаз между вращательными движениями соседних профилей, соответствующими крутильным колебаниям и равный при изгибных колебаниях - у2;
— сдвиг фаз между поступательным и вращательным движениями профиля - ф.
Отрицательным значениям у соответствовала волна деформации, бегущая по вращению колеса.
2.4. Результаты расчета устойчивости к флаттеру. Оценка влияния отдельный факторов
На рисунке 5 представлено изменение коэффициента аэродинамического демпфирования в зависимости от сдвига фаз между соседними профилями для изгибных, крутильных и изгибно-крутильных колебаний для двух режимов работы с расходами О1 и О 2.
Как видно из графика, коэффициент аэродинамического демпфирования отрицателен при крутильных колебаниях со сдвигом фаз у от -15° до -90°. Минимальное демпфирование было получено при крутильных колебаниях с межлопаточным сдвигом фаз = -60°. При этом расчеты с разными значениями расхода дали близкие результаты.
В остальных расчетах, как и ранее, наблюдается аэродинамическое демпфирование колебаний профилей.
Расположение кривых на рисунке 5 показывает, что по сравнению с результатами для изгибных колебаний, изгибно-крутильные меньше демпфируются потоком, а чисто крутильные дают отрицательное значение коэффициента аэродемпфирования, что создает условия для возникновения флаттера. Это значит, что, если при некоторых условиях колебаний лопаток в системе с диском, они имеют возможность совершать достаточно интенсивные крутильные колебания, то возможно возникновение флаттера.
Также в работе были проведены исследования условий возникновения флаттера в решетке профилей, совершающих изгибно-крутиль-ные колебания по первой форме с увеличенной относительно исходной в 1.3 раза амплитудой кручения. Такое значение представляется более обоснованным вследствие влияния диска (см. рис. 2).
Полученные зависимости коэффициента аэродемпфирования от сдвигов фаз между перемещением и вращением соседних профилей у при фиксированных сдвигах фаз между поступательным и вращательным движениями профиля Ф =-45° и ф = 0° представлены на рисунке 6 для двух режимов обтекания с приведенными осевыми скоростями в решетке и
Рис. 5. Изменение коэффициента аэродемпфирования в зависимости от сдвига фаз между колебаниями соседних профилей
Рис. 6. Зависимости коэффициента аэродемпфирования от сдвига фаз ф при фиксированных у = 0° и у = -30° для » = 0,00325
Анализ графиков показывает, что коэффициент аэродемпфирования принимает минимальные значения для всех расчетных случаев при у от -30° до 0°. При этом, в каждом расчете при уменьшении приведенной осевой скорости 1, коэффициент аэродемпфирования уменьшается.
На рисунке 7 изображены поля давлений в решетке профилей, совершающих изгибно-кру-
тильные колебания со сдвигом фаз между соседними профилями у = -30° за период Т. Стрелкой показано полученное направление движения волн. В эксперименте флаттер возникал при скорости с бегущей по вращению колеса волной деформации с двумя узловыми диаметрами, что соответствует у = -30°.
N = 0 ШШ1'
■р)
Щ 130000
■L щ
1
Л '
N = 23 В
t = 0-T
t = 0,2-T
t = 0,4-T
t = 0,6-T
t = 0,8-T
Рис. 7. Поля давлений в решетке за период колебаний Т
В работе также было проведено исследование влияния давления на входе в решетку на значение работы аэродинамического потока над лопаткой. Увеличение давления в 2 раза подразумевает также увеличение в 2 раза статической плотности при сохранении температуры на входе и выходе из решетки. Расчеты показали, что работа аэродинамического потока при такой постановке задачи увеличивается в 1,94 раза.
Заключение
1. С использованием разработанного метода и аэроупругих характеристик рабочих лопаток, характерных для первой ступени вентилятора, исследовано влияние отдельных факторов на устойчивость лопаток к флаттеру. Проведенные расчеты показали, что:
— при дросселировании ступени по напорной характеристике устойчивость лопаток к флаттеру снижается;
— при изгибных колебаниях лопаток всегда происходит их демпфирование потоком;
— при изгибно-крутильных колебаниях устойчивость к флаттеру снижается с увеличением крутильной составляющей колебаний;
— при учете влияния механической связанности лопаток через диск устойчивость к флаттеру снижается;
— существенное влияние на устойчивость лопаток к флаттеру имеет сдвиг фаз между лопатками в колесе и сдвиг фаз между изгибной и крутильной составляющими колебаний.
2. Валидация разработанного метода расчета проведена по имеющимся в ЦИАМ экспериментальным данным. При этом получено хорошее совпадение по границе флаттера, если эквивалентный профиль располагается на высоте 90% от максимальной длины лопатки, а расчет частот и форм колебаний производится с учетом влияния диска. При возникновении флаттера как в расчете, так и в эксперименте, имеет место волна деформации, бегущая по вращению колеса с двумя узловыми диаметрами, которая возникает при уменьшении расхода воздуха с постоянной частотой вращения ротора.
Литература
1. Хориков А.А. Прогнозирование и диагностика флаттера лопаток осевых компрессоров авиационных ГТД [Текст]/ А.А. Хориков // Труды ЦИАМ № 1311, Москва, 2012.
2. Программный комплекс Cobra v2.5. Свидетельство о государственной регистрации №2010613209 от 14 мая 2010г. Заявка №2010611632. Дата поступления 30 марта 2010 г.
Поступила в редакцию 01.06.2012
A.A. Khorikov, V.V. Shkurov, V.A. Shorstov. To the problem of compressor blisk blades flutter predicting
The calculation method of the blisk design impellers blades resistance to the bend-torsional flutter, based on the computation of the unsteady aerodynamic forces in the viscous flow for the oscillation period, has been proposed. The method reason is given on the basis of a mathematical flutter model, which takes into account the mechanical and aerodynamic relations of the blades in the impeller. To calculate the unsteady aerodynamic forces the developed in CIAM program complex Cobra was used. The validation of the proposed method is carried out on experimental data obtained from the fan test on stand C-3 CIAM. The consideration of the deformation disk influence, the phase difference influence between bending and torsional vibrations and the influence of the structural damping allowed to get good agreement with experiment.
Key words: blades, blisk, flutter-, frequencies and mode shapes, unsteady force, phase lag, energy method, damping, viscous flow.