CC BY
35
МЕХАНИКА
1
УДК 539.3+519.85+629.7.002
Ю. М. Темис, А. В. Троицкий
ПРОЕКТИРОВАНИЕ ОПТИМАЛЬНОГО ДИСКА ТУРБИНЫ
На примере задачи проектирования оптимального диска турбины показана эффективность применения метода анализа чувствительности. Рассмотрена задача минимизации массы диска при наличии прочностных и конструктивных ограничений. Метод позволяет находить оптимальные проекты дисков с учетом упруго-пластических деформаций. Приведен результат сравнительного анализа эффективности метода проекции градиента с использованием анализа чувствительности и метода квадратичной аппроксимации с применением конечноразностной аппроксимации производной. Показано, что метод анализа чувствительности позволяет существенно сократить время расчета и увеличить количество используемых параметров проектирования.
В современных газотурбинных двигателях применяются диски, конструкции которых создаются методами оптимального проектирования [1]. Оптимальным является диск минимальной массы, конструкция которого удовлетворяет заданным запасам прочности, технологическим и конструктивным ограничениям. При отсутствии градиентов температур, технологических и конструктивных ограничений возможно спрофилировать равнопрочный диск [1, 2]. Однако в общем случае создать равнопрочный диск, удовлетворяющий технологическим и конструктивным ограничениям, работающий при неравномерном распределении температуры по объему диска, невозможно. Поэтому в работах [1, 3] была развита методика профилирования оптимальных дисков, удовлетворяющих заданным запасам прочности, технологическим и конструктивным ограничениям. Применявшиеся в работе [1] методы решения задач оптимального проектирования относились к методам нулевого порядка. При этом время решения задач зависело от числа параметров проектирования, а качество проекта зависело от способности алгоритма приблизиться к глобальному минимуму при наличии ограничений. В работе [3] предпринята попытка на основе анализа моделей, описывающих предельное состояние диска, снять ограничение на число параметров проектирования. Однако алгоритм профилирования имел при этом существенные ограничения в возможности учета конструктивных ограничений и реального распределения напряжений и деформаций по объему диска.
В настоящее время разработаны алгоритмы оптимизации, позволяющие ослабить ограничения на число управляющих параметров в задаче. В основе этих алгоритмов лежит рассматриваемая в настоящей работе комбинация метода анализа чувствительности проекта [4, 5] к изменению управляющих параметров и методов минимизации первого или второго порядка точности.
Напряженно-деформированное состояние диска определяется следующей системой дифференциальных уравнений [1,2]:
¿и ¿т
а 1 Ти + —БТ
т Ет
а21
н
Я + а(1 + а)т,
= ЬНи + + (-ЕаТ - ри2т2)Н,
ат т т
(1)
где Я = аг Нт; т — радиус диска; аг — распределение радиальных напряжений по объему диска; и — перемещения точек диска; Е — модуль упругости материала диска; а — коэффициент Пуассона; Т(т) — температура точек диска; и — угловая частота вращения диска; р — плотность материала; а — коэффициент температурного расширения; Н(т) — толщины точек диска.
Решение системы (1) получим для краевых условий общего вида
ахи(а) + вх Я (а) = Пх,
а2 и(Ь) + в2 Я(Ь) = П2;
здесь ах, а2, в, в2 — некоторые коэффициенты.
Линейные соотношения (2) определяют способы задания граничных условий, приведенные в таблице. Заметим, что возможна комбинация любого столбца для т = а с любым из столбцов для т = Ь.
(2)
Задано перемещение
Задана сила
Задана упругая связь
а
1
А
-1/а
П1
0
= ь
Задано перемещение
Задана сила
Задана упругая связь
а
1
кь
в
П2
-1/Ь
иь
N Ь
г = а
0
к
п
0
1
и
а
Г
0
0
1
0
Целью оптимального проектирования является минимизация массы диска
^ [к] = 2п рк(г)гйг
(3)
при наличии ряда ограничений. Например, в некоторых точках интервала [а, Ь] могут быть наложены ограничения на толщины:
кшт(г) < к(г) < кшах(г). (4)
Ограничения на интенсивность напряжений представим в виде
(г) = ^О? + - Ог< [о]
где [о] — допускаемое напряжение в точке на радиусе
О О и ог = -—, = V---1— Е — Е аТ.
кг
кг
Таким образом,
(5)
ог(г) = Л 1(О)2 + (VО + иЕ — ЕаТ)2 — О- (VО + и-Е — ЕаТ\ . кг кг г кг кг г
(6)
Представим краевую задачу (1), (2) в операторном виде:
[з ]
{й — {й
= 0, [3] =
¿г
аг
—Ьг
а2
йг
Ь2
а1х1(а) + @гх2 (а) = щ, а?Хг(Ь) + в?Х2(Ь) = П?,
г £ [а, Ь],
— 1 — —2
где хг = и, х2 = О, аг =--, Ьг = ———, сг = а(1 + —)Т, а2
г Ег
Ь2 = —, с2 = (—ЕаТ — )к, = рш2г. г
Ограничение (5) заменим интегральным соотношением
(7)
Ек — '
««-/([ (ой2
1
.у
йг = Ф1(х1,х2 ,к,г)йг = 0, (8)
где [/]+ = тах{0, /} для любой функции /.
ь
г
г
ь
Оценим чувствительность функции цели (3), оператора задачи (7) и ограничения (8) к изменению параметра проекта, в качестве которого рассматривается толщина диска Н(г).
Рассмотрим два проекта дисков, имеющих толщины Нх (г) и Н2(г) = = Нх(г) + V (г).
Масса диска толщиной Н2 будет отличаться на 6^ [V] от массы диска толщиной Нь
SF [v ] = 2п prvdr = wvdr, w = 2 npr. (9)
Функционал Fi также изменяется при переходе от проекта h к проекту h1 + v, что в предположении малости v можно представить в виде
b
fid Ф1 д Ф1 д Ф1 \ SFi[v] = Fi[hi + v] - Fi [hi] ъ / i^Sxi + -^Sx^ + -dhv)dr, (10)
a
где {Sж}т = {Sxi, Sx:} — вариация вектора {ж}т = {x1 , x:}.
Непосредственно оценить SF1 при варьировании толщины диска не удается, так как в подынтегральное выражение явно входят вариации Sx1 переменных состояния xi, i = 1, 2, для определения которых необходимо решать краевую задачу (7) при измененном профиле hi (r) + + v(r). Отметим, что переход от вектора состояния {x} к вектору {x + + Sx} следующим образом изменит операторное уравнение (7):
л {:+а -й+--j {x:} Sh - ш Sh=о. с»)
Исключая из равенства (11) уравнение, соответствующее начальному состоянию ({ж}, Н), получим уравнение в вариациях и граничные условия:
(д/Л
j {S x:}+<
- h - f:
> Sh = 0, (12)
-h
a1Sx1(a) + eiSx: (a) = 0, a:Sxi(b) + P:Sx:(b) = 0,
b
b
где
ГдМ
д Н /
у д Н >
> =
А[7 ]! ± 1СЛ
д Н \х2\ д Н\с2\
(14)
Умножим выражение (12) на сопряженный вектор {ф}т = {фх, ф2} и проинтегрируем на отрезке [а, Ь]:
{ф}т
/ \д /1 ] \
д Н
7]{5х} + < > 5 Н
д /2
V . д Н , /
йт = 0.
(15)
Первое слагаемое в подынтегральном выражении можно проинтегрировать по частям, что позволяет представить равенство (15) в виде
Ь Ь _ _
д/- ВЦ
{5х}т [7 ]*{ф}йт +
где
[7 Г =
д Н' д Н
й йт -Ь-
{ф }5 Нйт + {ф}т {5х}
= 0, (16)
а-
-а2 Л Ь
---Ь2
йт
(17)
[7 ]'
ш
=
Введем сопряженную систему уравнений
д Фх д хх
д Ф-. д х 2 >
Тогда из равенства (16) следует, что соотношение (10) можно заме-
.
(18)
нить зависимостью
ОД [V] =
ь
/{ж
д /2 д Фх 1 ]
фх + "д^ф2 + ~дк\ V(т)йт = Мт»(т)йт, (19)
Г д/- / дФх \
"1(т) = | Ж, Ж, Ж/ у*
= 0. С учетом
при условии, что в выражении (16) имеем {ф}т{5х} условий (13) это налагает на решение системы (16) следующие гранич ные условия:
- ¡З-ф^а) + а-ф*(а) = 0,
- + а2ф2(Ь) = 0.
(20)
ь
а
ь
Для определения функции v необходимо решить задачу поиска минимума линейного функционала
ь ь
Е [Н х] + 6Е[v] = 2п / рНг(г)Г(1Г + w(г)v(Г)(1Г =
= C + J w(r)v(r)dr ^ min,
a
с интегральными ограничениями (21)
b
Fi[hi] + w1(r)v (r)dr = 0,
J (V(г))2dг — е2 = 0,
а
где С — масса исходного проекта диска; е — ограничения, налагаемые на норму изменения Н(г) на каждой итерации.
Задачу (21) с использованием множителей Лагранжа А, 7 можно свести к задаче поиска стационарной точки функционала:
b
b
ь ( ь
Ь[и] = С + [ т(г)v(г)dг — А ^ Ех + [ Wl(г)v(г)dг+
(22)
+2 (v(г))2dг — е2
Изменение функционала Ь^] при переходе от функции v к v + 6^] с точностью до линейных слагаемых выражается следующим образом:
ь ь ь
Ь[и + 6v] ^ Ь[и] +J w6vdг — А У тх6vdг + JJ v6vdг = Ь[и] + 6L[6v].
а а а
(23)
Как и в случае функции многих переменных, в стационарной точке линейная часть изменения функционала (вариация функционала) должна быть нулевой, т.е.
ь
6 Ь[6 v] = (т — Атх + Yv )6vdг = 0.
Из этого выражения следует
1
v = — (-w + Xwi).
Y
Для определения коэффициента А это выражение подставляется во второе из условий (21). Окончательное выражение для вариации толщины диска на каждой итерации имеет вид
/
v(r) = S
\
wwidr
—w +
Wi
wYdr
\
Fi
wi;
(25)
wYdr
/
здесь S = 1/y — параметр, который совместно с параметром г определяет норму изменения проекта на каждой итерации. Поиск оптимального проекта проведем методом проекции градиента, алгоритм которого представлен на рис. 1.
Изложенная методика применима только для задач без конструктивных ограничений на толщину диска. В случае наличия ограничений на минимальную и максимальную толщину диска произведем замену переменных
h =
km
+ k ■ k — k ■ " "'mm "'max ' bmin
--\--sin p
2 2 F
и решим задачу оптимизации относительно нового управления р.
Метод проекции градиента легко обобщается для случая наличия произвольного числа ограничений. При этом увеличивается число множителей Лагранжа, что несколько усложняет алгоритм их определения и выражение для расчета параметров изменения проекта.
В качестве примера рассмотрим случай двух прочностных ограничений: аг(г) < [аг] и а0(г) < [а0]. Как и в случае одного ограничения, эти локальные ограничения преобразуются в изопериметрические:
F [hi ] =
ük] ч +)
dr = Фl(xl,x2,h)dr = 0, (26)
F [h2] =
] 1 +)
dr = ^2(xl,x2,h)dr = 0.
(27)
b
b
b
b
b
2
b
b
2
Рис. 1. Алгоритм метода проекции градиента
Необходимо дважды выполнить анализ чувствительности для функционалов (25) и (26). Введем сопряженные уравнения
у ]'
д Ф„
д х1
д Фк
д Х2
к = 1, 2,
(28)
с граничными условиями (20), вариации функционалов запишем в виде
= 'ш1{г)у{г)в,г, ^2^] = w2(r)v (г)(г, (29)
д Ф1 д Ф1 т д Ф1 (1)
где = + —— 1р1' + —— 1ру2 '
дФ2 , дФ2 ,(1) , дФ2 ,(2) т2 = -тгг + —— Щ + т;— Ф2 •
д к дх
1
дх2
дк дх1 дх2
Искомая вариация управления имеет вид, аналогичный выраже
нию (24):
к-
V = — | — т
1
+ У^ \тк
к=1
)
(30)
В этой зависимости коэффициенты определяются подстановкой выражения (30) для v в линеаризованные изопериметрические ограничения (26), (27). Для определения коэффициентов Х^ получим систему двух линейных алгебраических уравнений:
2 ь ь
^^ Хк J гш1{т)гшк(г) (г = J гш1(г)гш(г)(г — ^Р1[К],
к=1 2
(31)
хк ! w2(r)wk(т)(1г = j w2(r)w(r)dr — jf2[к].
к=1
Рассмотренный алгоритм с помощью метода переменных параметров упругости можно обобщить для оптимизации диска, материал которого нагружен выше предела упругости, или диска из нелинейно упругого материала. При применении метода переменных параметров упругости рассматривается эквивалентное неоднородное упругое тело. Итерационная процедура расчета этим методом предполагает, что формально соотношение, описывающее связь между деформациями и напряжениями, остается таким же, как и в случае линейной задачи, но под модулем Юнга и коэффициентом Пуассона понимаются другие величины Е' и ^' [1], которые являются функциями координат и которые
ь
ь
ь
ь
можно найти в результате итерационного процесса решения физически нелинейной задачи. Интенсивность напряжений является функцией интенсивности деформаций и определяется обобщенной диаграммой деформирования а = ), е^ — интенсивность деформаций.
Учет пластичности методом переменых параметров упругости приводит к необходимости на каждой итерации задачи оптимизации итерационным методом определять значения вектора {ж} в точках диска. Полученные в результате решения нелинейной задачи характеристики материала Е и ^ используются при решении сопряженной задачи.
Для решения исходной краевой задачи (7) и сопряженной системы (18), (20) удобно использовать метод начальных параметров. В этом случае для исходной задачи интегрируется система дифференциальных уравнений (7) с нулевым столбцом свободных членов {с} при начальных условиях {ж(1)(а)}т = {0 1} и {ж(2)(а)}т = {1 0} для получения двух линейно независимых частных решений однородного уравнения. Далее интегрируется система (12) с неоднородной правой частью (ненулевым столбцом {с}) для получения частного неоднородного решения {ж(0)(г)} с произвольными начальными условиями. Эти три решения позволяют получить общее решение задачи (12), (13) в виде
Для получения констант Ск соотношение (31) подставляется в граничные условия задачи (13). Сопряженная задача решается аналогично.
Эффективность применяемого метода анализа чувствительности оценивалась по результатам сравнения его с методом квадратичной аппроксимации (методом спроектированных лагранжианов) [6, 7], использующим традиционное численное дифференцирование целевой функции и ограничений. На рис. 2 представлен алгоритм метода квадратичной аппроксимации [6], применяемый при решении общей задачи оптимизации для п параметров проектирования Л, целевой функции f (Л) и т нелинейных ограничений дз(Н)\
2
(32)
к=1
тах
г = 1, 2, ... ,п;
з = 1, 2,...,те;
к = те + 1,... ,т.
(33)
дз(Л) = ° дк (Л) < 0,
Конец
Рис. 2. Алгоритм оптимизации методом квадратичной аппроксимации
Общая схема системы оптимального проектировании конструкции диска представлена на рис. 3. На этой схеме стрелками обозначено направление потока данных между отдельными блоками. В блоке анализа чувствительности проекта в общем случае для расчета градиентов мож-но применить численное дифференцирование, проводя серию расчетов напряженно-деформированного состояния для вычисления неявно заданных через выходные параметры целевой функции и ограничений. При использовании более совершенных методик анализ чувствительности тесно связан с расчетом напряженно-деформированного состояния. В связи с этим в алгоритме оптимизации производные вычисляются непосредственным численным дифференцированием, либо, если имеется такая возможность, их значения берутся из блока прочностного расчета конструкции. При профилировании диска методом квадратичной аппроксимации для расчета градиентов целевой функции и ограничений использовалось непосредственное численное дифференцирование. При использовании такой методики проводится п расчетов напряженно-деформированного состояния для получения всех необходимых производных, что приводит к возрастанию продолжительности расчета и, соответственно, накладывает существенное ограничение на число управляющих параметров.
Рис. 3. Общая схема системы оптимального проектирования (здесь
д (к) = {91 (к),...,д, (к)})
Сравнение двух методов проводилось при решении задачи оптимизации диска с учетом линейного изменения температуры (550 0С у ступицы и 650 0С у обода), ограничений на максимальную (70 мм) и минимальную (5,4 мм) толщину диска, переменных допускаемых напряжений (1050 МПа у ступицы и 740 МПа у обода) с внутренним радиусом а = 108 мм, внешним Ь = 288 мм и кь = 46 мм; частота вращения составляла п = 10700 об./мин.
На рис. 4 представлен оптимальный проект диска, полученный методом проекции градиента с 200 управляющими параметрами, на рис. 5 — результат оптимизации диска при тех же нагрузках с дополнительными конструктивными ограничениями (штриховыми линиями показано начальное приближение).
Рис.4. Оптимизация диска методом проекции градиента: 1 — аг, 2 — а^,
3 — а^, 4 — [а]
Рис. 5. Оптимизация диска с дополнительным конструктивным ограничением методом проекции градиента:
1 — аг, 2 — а^, 3 — , 4
[а ]
Ь,ММ 60 4/7 20 О
108
1 1 —
/ К
1
т 180
116 252 г, мм
Рис. 6. Оптимальный проект диска , полученный методами проекции градиента (1) и квадратичной аппроксимации (2)
Рис. 7. Оптимальный проект диска с дополнительным конструктивным ограничением, полученный методами проекции градиента (1) и квадратичной аппроксимации (2)
Решение задачи профилирования диска методом проекции гради-ента (25) и методом квадратичной аппроксимации показало, что результаты расчетов незначительно отличаются как по массе, так и по форме конструкции (рис. 6, 7). Однако следует заметить, что применение полуаналитического анализа чувствительности проекта позволило значительно сократить продолжительность расчета. При применении метода квадратичной аппроксимации использовалась параметризация диска в виде ломаной из 30 отрезков. Это привело к необходимости проводить расчет конструкции 30 раз для вычисления градиентов на каждой итерации. Таким образом, продолжительность расчета в этом случае существенно зависела от числа управляющих параметров. В методе проекции градиента с полуаналитическим анализом чувстви-тельности на каждой итерации проводится два расчета напряженно -
деформированного состояния, независимо от числа управляющих параметров, при этом расчет градиентов производится без погрешностей, присущих методам численного дифференцирования. Таким образом, очевидно, что продолжительность расчета при применении метода квадратичной аппроксимации без анализа чувствительности приблизительно в 15 и более раз превосходит продолжительность расчета при применении метода проекции градиента в данной постановке задачи. Для параметризации диска при использовании метода проекции градиента использовалась ломаная из 200 и более отрезков, при этом продолжительность расчета была значительно меньше, чем при использовании квадратичной аппроксимации с 30 параметрами.
Из графиков изменения массы диска по итерациям, приведенных на рис. 8, видно, что метод проекции градиента при решении данной задачи имеет более высокую скорость сходимости по сравнению с методом квадратичной аппроксимации. Это объяснятся тем, что в реализованном методе проекции градиента величина шага выбиралась пользователем и в дальнейшем уменьшалась при необходимости. В методе квадратичной аппроксимации величина шага выбиралась из условия минимума функции выигрыша [6], что обеспечивает более надежную сходимость к истинному решению, но может приводить к вычислению малых шагов в пространстве переменных проектирования.
Таким образом, полуаналитический метод анализа чувствительности позволяет существенно увеличить эффективность алгоритма опти-
10 65
\ 60 ъ
(о <о
§ 55 50 ¿>5
О 5 10 15 20 30 ио 50 60 80 100 120 140 160 Номер итерации
Рис. 8. Зависимость массы диска от итерации при использовании методов квадратичной аппроксимации (1) и проекции градиента с различными значениями начального шага (2-4) (значение шага максимально для кривой 4)
мизации, требующего расчета производных. Значительное уменьшение продолжительности расчета позволяет использовать большое количество параметров проектирования и получать более точные решения задач оптимизации.
СПИСОК ЛИТЕРАТУРЫ
1. Демьянушко И. В., Биргер И. А. Расчет на прочность вращающихся дисков. - М.: Машиностроение, 1978. - 247 с.
2. Демьянушко И. В., Королева Е. Ф. Оптимальное проектирование дисков турбомашин // Изв. АН СССР. Сер. Механика твердого тела. - 1972 - № 2. - С.176-180.
3. Братчик В. Я. Численное проектирование дисков ГТД методом динамического программирования // Труды ЦИАМ № 996: Проблемы прочности и динамики в авиадвигателестроении. - 1982. - Вып. 2.
4. Х о г Э., А р о р а Я. Прикладное оптимальное проектирование. - М.: Мир, 1983.-479 с.
5. ХогЭ., Чой К., Комков В. Анализ чувствительности при проектировании конструкций. - М.: Мир, 1988. -428 с.
6. Г и л л Ф., МюррейУ, Райт М. Практическая оптимизация. - М.: Мир, 1985.-509 с.
7. F l e u r y C. First and second order convex approximation strategies in structural optimization // Structural Optimization. - 1989. - № 1. - P. 3-10.
Статья поступила в редакцию 18.02.2003
Юрий Моисеевич Темис родился в 1947 г., окончил в 1971 г. МВТУ им. Н.Э. Баумана. Д-р техн. наук, профессор кафедры "Прикладная математика" МГТУ им. Н.Э. Баумана, руководитель отдела математического моделирования Центрального института авиационного моторостроения им. П.И. Баранова (ЦИАМ). Автор более 100 научных работ в области прочности и динамики турбомашин, численных методов решения нелинейных задач упругости, пластичности и ползучести.
Yu.M. Temis (b. 1947) graduated from the Bauman Moscow Higher Technical School in 1971. D. Sc. (Eng), professor of "Applied Mathematics" department of the Bauman Moscow State Technical University and chief of mathematical simulation department of the Central Institute for Aviation Motor Engineering (TslAM) n. a. P.I. Baranov. Author of over 100 publications in the field of machine strength and dynamics, numerical methods of solving nonlinear problems of elasticity, plasticity and creepage.
Артем Владимирович Троицкий родился в 1980 г., окончил в 2003 г. МГТУ им. Н.Э. Баумана. Аспирант кафедры "Прикладная математика" МГТУ им. Н.Э. Баумана.
A.V. Troitsky (b. 1980) graduated from the Bauman Moscow State Technical University in 2003. Post-graduate of "Applied Mathematics" department of the Bauman Moscow State Technical University.
