Ресурсная оптимизация деталей газотурбинного двигателя Текст научной статьи по специальности «Механика и машиностроение»

УДК 534.1:539.3

РЕСУРСНАЯ ОПТИМИЗАЦИЯ ДЕТАЛЕЙ ГАЗОТУРБИННОГО ДВИГАТЕЛЯ

1 2 Нгуен Динь Дыонг1, В.И.Рыжиков2

Иркутский государственный технический университет,

664074, г. Иркутск, ул. Лермонтова, 83.

Статья содержит краткое описание методов для динамики и прочности лопаток турбомашин, расчетов чувствительности и оптимального проектирования лопаток. Результаты исследований позволяют сделать вывод о том, что анализ чувствительности представляет собой весьма эффективное средство при проектировании сложных машиностроительных конструкций. Ил. 11. Табл. 3. Библиогр. 9 назв.

Ключевые слова: ресурс деталей турбомашин; метод конечных элементов; оптимизация проектирования турбомашин; чувствительность динамических напряжений; долговечность.

RESOURCE OPTIMIZATION OF GAS-TURBINE ENGINE PARTS Nguyen Dinh Duong, V.I. Ryzhikov

Irkutsk State Technical University, 83 Lermontov St., Irkutsk, 664074.

This article contains a brief description of methods for the dynamics and the strength of turbomachine blades, calculations of sensitivity and optimal design of turbine blades. Testing results suggest that the sensitivity analysis is a very effective tool in designing complex engineering structures. 11 figures. 3 tables. 9 sources.

Key words: resource of turbomachinery parts; finite element method; optimization of turbomachinery design; dynamic stress response; durability.

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

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

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

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

Основные математические модели описываются уравнением движения Лагранжа. В общем случае на-гружения конструкции имеют вид [2]:

МЁ + СЁ + КЁ = ^ ((), (1) где - вектор внешних сил; С - матрица демпфирования; М - матрица масс конструкции; К - матрица жесткости конструкции; 1 - вектор перемещений конструкции.

При работе двигателя сила, действующая на всю структуру, может быть получена как сумма конечного числа гармоник Фурье для каждой подструктуры [3]. Сила Рк] действует В направлении к-й степени свободы ]-й подструктуры конечно-элементной модели рабочего колеса в момент времени t и записывается в виде конечного ряда Фурье [1, 3]:

N

Рщ ) = Е Арк () • е« (2)

р=1

где ЛркА) - текущее значение р-й рассматриваемой гармоники Фурье от сил, соответствующих константе у волнового распространения; N - число подструктур.

1Нгуен Динь Дыонг, аспирант, тел.: 89246038970, e-mail: nguyendihnduong1983@yahoo.com Nguyen Dinh Duong, postgraduate student, tel.: 89246038970, e-mail: nguyendihnduong1983@yahoo.com

2Рыжиков Вячеслав Игоревич, аспирант, тел.: 89086629213, e-mail: x_mybox@mail.ru Ryzhikov Vyacheslav Igorevich, postgraduate student, tel.: 89086629213, e-mail: x_mybox@mail.ru

Уравнение (2) соответствует распределению сил, действующих на структуру. Из (2) имеем [4]:

Apk (t) = (1/ N) -X Fkj (t) • е

- i •( - j )-м

(3)

Каждая из рассматриваемых гармоник Фурье приобретает в соответствии с уравнением (3) изменение на каждой лопатке. В итоге уравнение силовой вибрации у-й подструктуры для конкретных гармоник Фурье для переходных режимов имеет вид

Мр -¿З (I) + Ср -¿р (I) + Кр ■ ¿¿> (I) = Ар (О, (4) где р=1,2,..,М, а Мр, Ср и Кр - матрицы масс, демпфирования и жесткости подструктуры для р-й гармоники соответствующих сил и & - вектор перемещений подструктуры. Вектор А (V) включает все значения

Арк (V).

Для решения уравнения (4) для к-й подструктуры с изменением времени I получаем: для N - четного

dk (t) = 3J (t) + 2:

N/2 p=1

(5)

3N/2(t)+2 • X

p=1

(6)

Re (p (t)) • Cos(k - J) • (2np / N) +

+Im(3(t)) • Sin(k - J) • (2np / N) для N - нечетного dk (t) = 33 (t)+(-1)k-1 •

Re (33 (t)) • Cos(k - J) • (2np / N)

+Im(33 (t)) • Sin(k - J) • (2np / N)

где dk(t) - общий вектор перемещений k-й подструктуры в момент времени t.

Из этих уравнений общее решение для силовой вибрации на переходных режимах получается из выражений (5) и (6) при соответствующих значениях р = 0,1,2,..., (N-1)/2 для N - нечетного, и р = 0,1,2,...N/2 для N - четного.

Для требуемого числа групп m уравнение движения примет вид

Pi + 2<1 'Щ -Pt + ®2iPi = fi(t), (7)

гдеpl = pl (t);щ =щ (t), l=1,2,..m; < -модальный

коэффициент демпфирования.

Для силы f0, действующей в течение периода üt, отклик системы (модальные перемещения) имеет вид

[5]:

fC1 -Cosa,- V(1 -<)• At + А

+C2 • Sinat • ^(1 -<) • At + (Fо /af)y

p =

где C1 =

C 2 =

XQ -(F0 / u? )

V0 • ®I •(x0 -(F0 /а 2 )

t- V(1-<t2)

а Х0 и У0 - модальные перемещения и скорость в начале периода Д1, модальный коэффициент демпфирования.

Модальные силы возбуждения 1 имеют вид [6]:

1 (I) = ЗТрГАр (I).

Реальные перемещения получаем из выражения

[6]:

т

¿(V) ^¿Т'Р, (0, 1=1,2,т

1=1

Модальный отклик подструктуры (перемещения, напряжения) определяется из уравнения (7). Общий отклик полной системы определяется в соответствии с уравнениями (5) и (6).

При решении задачи оптимизации в [2] целевой функцией является первая частота собственных колебаний, а варьируемыми параметрами - толщины профиля лопатки в узлах конечноэлементной сетки. Собственные частоты и формы колебаний определяются при решении задачи

(К -а2М)u = 0,

(8)

где К - матрица жесткости; М - матрица масс; ш - собственная частота колебаний; и - собственная форма колебаний, нормированная условием

(Ми, и) = 1. (9)

Для решения проблемы собственных значений (8, 9) используется алгоритм Ланцоша.

Вычисление градиента целевой функции ш. При изменении толщины Ьк в каком-либо узле номер к изменяются матрицы жесткости и масс конечных элементов, прилегающих к этому узлу. Это влечет изменение общих матриц К и М и как следствие изменение собственных частот и форм. Продифференцируем соотношение (8):

дК 2 дМ _ да ля --а2--2а-M

dhk dhk dhk _

du

и + (K -а2 M)— = 0.

дК

Обозначим особенную матрицу динамической жесткости, вычисленную на собственной частоте:

2 = (К -а2М). (10)

Умножим выражение (10) скалярно на и. Второе слагаемое в силу симметричности матрицы динамической жесткости 2 и в силу (8) обращается в нуль:

, ди

-,и

fди „ ^

-, Zu

дК

= 0,

поэтому в (10) остается только первое слагаемое:

ff К

- - а

дМ

и,и

к /

_ да ... . _

- 2а-(Ми,и) = 0.

дК

Учитывая условие нормировки (9), получим формулу для производной от собственной частоты колебаний по толщине:

да

Ж

J_ 2а

fi

дК 2 дМ -а

дК

дК

и, и

к

= 0.

дК / дИк, дМ / дИк нахо-

Частные производные

дятся, как и сами матрицы К и М, путем суммирования соответствующих матриц отдельных элементов. При

к=1

этом для дК /дкк, дМ/дкк вычислительный процесс существенно ускоряется, так как нужно суммировать матрицы, производные только для тех элементов, которые прилегают к к-му узлу. Выражения для частных производных от К и М отдельного элемента получены в аналитическом виде путем дифференцирования выражений для элементов самих матриц К и М.

Методы оптимизации. Для решения задачи оптимизации в [9] используется метод, учитывающий информацию о градиенте. По мнению многих исследователей, наилучшие результаты дает применение квазиньютоновского метода (метода переменной метрики) с использованием формулы Бройдена-Флетчера Гольдфарба-Шенно (БРСБ-формулы) для пересчета гессиана [9].

Различные реализации этого метода отличаются, в

у2 и у3 выбирается минимальное - это и будет результат итерации.

уЛ

УI

х

-VI Л'з -?2

Рис. 1. Одномерный поиск

Если минимальным значением оказывается у1, то делается попытка возврата к диагональному гессиану, а если это не помогает, шаг I уменьшается вдвое.

Таблица 1

Форма колебаний п = 0 1/с п = 100 1/с

БТ!0 18 К1!БЛ [7] БТЮ 18 1Ч!БЛ [7]

1И 132.79 121 224.41 204

2И 455.15 421 545.09 528

1К 941.53 804 961.07 815

3И 1142.60 1048 1226.16 1138

2К 1638.34 1430 1697.39 1507

частности, разными схемами проведения одномерного поиска в найденном направлении. Если целевая функция не очень сильно отличается от квадратичной, то хорошие результаты дает классическая схема, когда проводится полный квазиньютоновский шаг:

б = - Н'1д, (11)

где д - градиент, а Н - аппроксимация гессиана на данной итерации. Однако при расчетах элементов машиностроительных конструкций целевая функция, как правило, весьма далека от квадратичной. Поэтому для одномерного поиска вдоль направления б можно использовать схему кубической интерполяции (рис. 1).

Задается начальная длина шага I, которая в процессе оптимизационного поиска будет уменьшаться или увеличиваться. Первая точка х1 - это результат предыдущей итерации. В этой точке известны значение целевой функции у1 и производная д1 (проекция градиента на направление б ). Если оказывается, что д1 > 0, то данное направление поиска признается неудачным. В этом случае надо вернуться к единичному гессиану и следующий шаг провести в направлении антиградиента. Там такая ситуация точно не возникнет.

Далее осуществляется шаг длиной I в направлении б и получается точка х2. В этой точке вычисляются функция у2 и градиент - находится д2. Если у2<у1 и при этом д2 < 0, то продолжается поиск в направлении б увеличивающимися шагами (I увеличивается в 2 раза), пока хотя бы одно из этих условий не будет нарушено. Теперь точно известно, что между двумя последними точками находится минимум кубической параболы -точка х3. В этой точке определяется целевая функция у3, ее градиент, и из трех значений целевой функции у1,

Таким образом, в процессе оптимизационного поиска длина шага I может как увеличиваться, так и уменьшаться, адаптируясь к поведению целевой функции. Критерий выхода - уменьшение шага I до заданной малой величины.

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

х~ < хк < х+,к е[1,п]. (12)

Этот метод основан на ограничении шага и скольжении вдоль ограничения. Для начала работы требуется стартовать из допустимой точки, которая может находиться внутри или на границе области, заданной ограничениями (12). С помощью уравнения (11) определяется направление поиска б и производится одномерный поиск, как описано ранее, но только до ближайшей границы (рис. 2,а).

Рис. 2. Учет ограничений

Если при этом точка х1 находится на границе, а вектор б направлен вдоль границы или внутрь области,

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

Численные результаты анализа чувствительности перемещений при колебаниях лопаток турбин ГТД. Для расчета реальных конструкций рассмотрена рабочая лопатка газовой турбины. Общий вид и

характеристики лопатки приведены на рис. 3. На рис. 4 показана конечно-элементная модель лопатки на основе треугольных оболочечных элементов с 6-ю степенями свободы (55 узлов, 80 элементов, 330 степеней свободы). Лопатка жестко защемлена в корне. Первые 4 формы колебаний лопатки изображены на рис. 5. Резонансная диаграмма лопатки в диапазоне 0...100 1/с для первых трех изгибных форм колебаний представлена на рис.6. Основные характеристики лопатки: = 0,3 м; I = 0228 м; а = 85,5°; у = 46,5°; Е = 2,1 х 105 МПа; р = 7,85 кг/м3; ц = 0,3. В табл. 1 приведены результаты расчета собственных частот лопатки [6].

Некоторое расхождение расчетов с помощью конечного элемента БТЮ 18 из пакета В1_А0!8+ и расчета на основе трехмерного конечного элемента из па-

55 44 33 22 11

52.

К

/\

7\

Z$

4WW

R

Рис. 3. Общий вид лопатки

/ (Гц) ■ 1000 1,6 1,4 1,2

1,0 0,8 0,6 0,4 0,2 0

Рис. 6. Резонансная диаграмма рабочей лопатки

28/

15/

14/

Рис. 4. Конечно-элементная модель лопатки

Сумма гармоник для одного периода

Рис. 5. Формы колебаний лопатки

Первая и вторая гармоники для одного периода

.0 .1 .2 .3 .4 .5

.9 1.0 1.1 1.2 1.3 1.4 1.5

Время (с) * . 1Е+00

-.2 -.4 -.6 -.8 -1.0 -1.2

.2 .3 .4 .5 .6 .7 .8 .9 1.0 1.1 1.2 1.3 1.4 1.5

Время (с) * . 1Е+0С

Рис. Т. Распределение 1-ой и 2-ой

гармоник возбуждения

W . .......

ГГШ :0 "W _

Рис. 8. Распределение суммарной нагрузки

Напряжения, Н/мм 2 .1Е+03 ( kt= 0,9)

1.0 2.0 3.0 4.0 5.0 6.0 7.0

Рис. 9. Динамические напряжения при демпфировании high

Рис. 10. Динамические напряжения при увеличении толщины

Рис. 11. Динамические напряжения при уменьшении толщины

26.0

25.0

24.0

8.0

кета КНБЛ [7] обусловлено различной степенью конечно-элементной идеализации лопатки и моделированием хвостовика в последнем случае.

Следующим этапом исследования явился расчет и анализ изменения собственных частот колебаний в зависимости от геометрических и эксплуатационных параметров. Влияние варьирования толщины всего профиля лопатки на частоты колебаний иллюстрирует табл. 2 к = 1.1 означает увеличение толщины всего профиля на 10%, а кл = 0.9 - уменьшение толщина на 10%).

Влияние изменения толщины всей.

Величина демпфирования для первых трех изгиб-ных форм и двух вариантов демпфирования - "low" и " high" - приведена в табл. 3.

Максимальные значения динамического отклика, а именно изгибные напряжения, были рассчитаны в центре тяжести треугольника с узлами 1,2,12 (см. рис. 4) при демпфировании "high". Исследовался разгон от 0 до 100 1/с в течение 25 с [7].

Рис. 9 показывает напряжения при демпфировании " high". Все дальнейшие расчеты были проведены для демпфирования " high". Рис. 10 и 11 иллюстриру-

Таблица2

Форма Кн = 1.1 Кн = 0.9

N= 0 1/c N= 100 1/c N= 0 1/c N= 100 1/c

f Af, % f Af,% f Af, % f Af, %

1И 142.62 +7.40 230.65 +2.78 122.75 -7.56 218.43 -2.66

2И 488.64 +7.36 572.36 +5.00 420.43 -7.63 517.34 -5.09

1К 992.74 +5.44 1010.69 +5.16 891.63 -5.44 913.21 -4.98

3И 1210.83 +5.97 1283.74 +4.70 1067.55 -5.97 1162.20 -5.22

2К 1695.27 +3.47 1756.00 +3.45 1580.28 -3.47 1639.98 -3.38

Уменьшение толщины лопатки ведет к уменьшению собственных частот и наоборот. Более сильно изменяются изгибные собственные формы, менее -крутильные.

В дальнейшем были исследованы вынужденные колебания для всех вышеописанных случаев. Лопатка возбуждалась через 20 сопловых лопаток (z =20). Фурье-ряд, например, для подъемной силы имеет вид [8]:

f (t) = L0(1 + 0,5 cos q> + 0,025 cos 2ф) , (13)

где первый член представляет собой статическую часть.

Вид этого возбуждения для одного периода показан на рис. 7, который изображает нагрузку от 2-го и 3-го членов уравнения (13), а на рис. 8 показана сумма статической и динамической составляющей нагрузки из уравнения (13). Силы возбуждения равномерно распределены по длине лопатки [6].

ют динамические напряжения при изменении толщины лопатки.

Таблица 3 Коэффициенты демпфирования [6]

Режим Значения коэффициентов

Z1 Z2 Z3

low 0.00075 0.00094 0.0011

high 0.00150 0.00190 0.0023

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

Библиографический список

1. Репецкий O.B. Компьютерный анализ динамики и прочности турбомашин. Иркутск, 1999. 301 с.

2. Заинчковский К.С. Разработка расчетно-оптимизацион-ных моделей для анализа прочностных и вибрационных характеристик лопаток ГТД: дис. ... канд. техн. наук. Иркутск, 1995. 136 с.

3. Omprakash V., Ramamurti V. Natural Frequencies of Bladed Disks by a combined Cyclic Symmetrie and Rayleigh-Ritz method //J. Sound and Vibration. 1988. V.25, N 2. P. 357366.

4. Downing S.D., Socle D.F. Simple rainflow counting algorithms // Int. J. Fatigue, January. 1982.

5. Ramamurti V. Computer Aided Design in Mechanical Engi-neering.-New Delhir Tata McGraw - Hill Publishers. 1989. 331 P.

6. Irretier H., Repetski 0., Sainchkowski K. Zu Empfindlich-keitsanalyeen, instationaren Schwingungen und Lebensdauern-

berechnungen von Turbinenschaufeln und Laufradern bei transienten und Stationaren Betriebsbedingungen. -Kassel:Kassel Universitat, Institut fur Meohanik, 1994. 26 S.

7. Hohlrieder M. Zur statischen und dynamischen Analyse rotierender elastischer Structuren (Turbinenschaufeln, Verdichter) bei transienten Betriebsbedingungen/ Dis. Dr. - Ing. Kassel Universitat: Institut fur Mechanik. 1994. 201 S.

8. Irretier H., Repetski 0. Analyse der Eigenschwingungen rotierender axialer und radialer Laufrader und Schaufelpakete von Turbomaschinen mittels Hyperelemente, Kondensation und der Methode zyklischer Symmetrie. - Vorschungsbericht DAAD, Mittellung 4/1991.- Kassel Universitat: GHK, Institut fur Mechanik, 1991. 121 S.

9. Оптимизация лопатки рабочего колеса турбокомпрессора по критерию собственной частоты колебаний / В.А.Жовдак [и др.] // Вестник. Харьков, 2003. № 12[1].