CC BY
23
УДК 621.438 - 253:531.3.001.57
Р.Л. ЗЕЛЕНСКИЙ, С.В. ЕПИФАНОВ, Е.В. МАРЦЕНЮК, В.В. БОЙКО
Национальный аэрокосмический университет им. Н.Е. Жуковского «ХАИ», Украина
МОДЕЛИРОВАНИЕ ДИНАМИКИ РАДИАЛЬНОГО РАЗМЕРА ДИСКА ТУРБИНЫ, ОБУСЛОВЛЕННОЙ ДЕЙСТВИЕМ МЕХАНИЧЕСКИХ ФАКТОРОВ
Важным параметром, влияющим на КПД турбины, является радиальный зазор между ротором и статором. Развитие систем активного управления зазором обусловливает необходимость разработки и использования математических моделей, представляющих динамику изменения зазора в различных условиях работы двигателя. Значение радиального зазора определяется радиальными размерами деталей, участвующих в его формировании: диска, лопатки и корпуса. В данной работе рассмотрено моделирование радиального размера диска. Для этого необходимо учитывать не только нагрев, но и действие силовых факторов нагружения. Использование упрощенных моделей соответствующих деформаций не дает приемлемой точности, так как упрощенная модель не учитывает изменение механических свойств (в частности, модуля упругости материала) от температуры. Предложена методика расчета радиального удлинения диска от действия силовых факторов с динамической коррекцией значения модуля упругости по температуре. Сравнительные результаты моделирования удлинения с использованием МКЭ и упрощенной модели с коррекцией модуля упругости показали, что предложенная модель обеспечивает необходимую точность.
Ключевые слова: газовая турбина, диск, динамическая математическая модель, радиальное удлинение, механическое нагружение, температура, переходные характери-
стики, граничные условия.
Введение
Развитие газотурбинных двигателей (ГТД) связано с постоянным интенсивным поиском резервов дальнейшего повышения эффективности рабочего процесса с целью улучшения их экономичности. Одним из таких резервов является сокращение потерь, связанных с перетеканием рабочего тела в различных зазорах между вращающимися и неподвижными деталями двигателя.
Опыт эксплуатации, а также многочисленные экспериментальные исследования убедительно показывают, что размеры радиальных зазоров (РЗ) между статором и ротором существенно влияют на КПД турбины и, следовательно, на удельный расход топлива ГТД. Установлено [1], что увеличение радиального зазора по торцам рабочих лопаток турбины на 1% вызывает уменьшение ее КПД, соответствующее возрастанию удельного расхода
топлива Суд на 1,52%. Еще большее влияние РЗ наблюдается в компрессоре — при его увеличении на 1% происходит возрастание на 1,53%.
К основным факторам, определяющим величину РЗ, относятся силовые воздействия на статор и ротор двигателя, а также термиче-
ское расширение их элементов на переходных режимах. Данный процесс будем называть прогревом двигателя.
Создание и эксплуатация сложных технических систем на современном уровне предполагает обязательное применение их математических моделей, отражающих информацию о системе [2].
Современный ГТД представляет собой сложную нелинейную динамическую систему с взаимным влиянием газодинамических и те-плофизических процессов, протекающих в его узлах. Процессы в двигателе являются нестационарными по времени и условиям эксплуатации, а для отдельных конструктивных схем имеют переменную структуру. Функционирование двигателя происходит при постоянном действии внутренних и внешних возмущений, а для маневренных самолетов — практически на переходных режимах.
При создании и эксплуатации ГТД разрабатывается и используется значительное количество математических моделей, описывающих свойства как двигателя в целом, так и его отдельных узлов. Это модели напряженно-деформированного, теплового состояния лопаток, дисков, роторов и других элементов компрессора и турбины, камеры сгорания, сопла, термогазодинамические модели, описывающие
© Р.Л. Зеленский, С.В. Епифанов, Е.В. Марценюк, В.В. Бойко, 2015 - 80 -
рабочий процесс в элементах двигателя, т.е. связи между давлением, температурой, расходом воздуха и газа в различных сечениях проточной части двигателя и другие модели.
В зависимости от решаемых задач управления и динамики ГТД применяются различные математические модели, оптимальные по степени адекватности описания процессов в двигателе, простоте и наглядности. В них учитываются различные факторы, влияющие на характеристики ГТД как объекта управления и его динамику; они различаются математической записью уравнений и методами их численного решения.
В данной работе рассмотрена задача формирования динамической математической модели рабочего процесса двигателя, учитывающей изменение степени прогрева конструктивных элементов.
Проблема моделирования динамики двигателя рассматривается в многочисленных статьях и монографиях [1 - 7]. Для моделирования процесса прогрева двигателя поузловую динамическую модель необходимо дополнить моделью влияния его прогрева на изменение параметров рабочего процесса. Прогрев, происходящий вследствие теплообмена между горячим газом и элементами конструкции, сопровождается изменением радиальных зазоров в лопаточных машинах, приводящим к изменению их характеристик.
Для определения зазоров необходимо разработать быстросчетную модель деформированного состояния ротора и статора двигателя, учитывающую начальную температуру деталей и динамику прогрева. Такая модель может быть выполнена на основе методических подходов, предложенных в работе [8]. Модель зазора реализована как разность радиальных перемещений ротора и статора. Перемещения ротора и статора турбины определяются из расчетов соответствующих деталей и их сборок при помощи пакета ANSYS.
1. Математическая модель перемещений
Рассмотрим узел конструкции двигателя, состоящий из неравномерно нагретых деталей,
нагруженных поверхностными силами P(x,y,z) и массовыми силами F(x,y,z). Созданные этими факторами перемещения u(x,y,z), деформации s(x,y,z) и напряжения a(x,y,z) в точке
(x,y,z) должны удовлетворять уравнениям теории упругости [9]
e = Ru ; (1)
Ятсг + Р = 0 (3)
с соблюдением граничных условий
Р-С<г = 0, (4)
где е = [ех Бу е2 Уху Уу.,. Угх ]т - вектор-столбец
деформаций; <у = [стх <уу ах сг^, а^ а^ ]т —
вектор-столбец напряжений; и = [ихиуи2]т — вектор перемещений;
"а/Эх 0 0 д/ду 0 д/дх
Я= 0 д/ду 0 Э/Эх д/дг О
_ 0 0 д/дт. О д/ду д/дх_
матричный дифференциальный оператор; ~1-ц ц ц О О О ц 1-ц ц О О О ц ц 1-ц 0 0 0 —матрица 1-ц-2ц* ООО 0.5-ц 0 0 0 0 0 0 0.5-ц 0 0 0 0 0 0 0.5-ц
модулей упругости для изотропной среды; Е — модуль упругости, ц — коэффициент Пуассона; а = [ах ау а2 аху а^ аж ]т — вектор коэффициентов температурного расширения; С — матрица направляющих косинусов пу 0 0 пу о п, 0 п, I
D = -
С =
0 п„
*z
о
o = D(s-at);
(2)
В соответствии с подходом, принятым в МКЭ, разобьем рассматриваемый узел конструкции на п трехмерных четырехточечных симплекс-элементов (тетраэдров). Применим наиболее разработанный и распространенный вариант МКЭ — метод перемещений [9], который представляет зависимость перемещения
и(е> в пределах элемента от координат в виде линейной аппроксимирующей функции
и(х, у, г) = N1 (х, у, г)и1 + N4 (х, у, г)и4 = N(e)u(e),
где — (3 х 3)-базисная функ-
ция аппроксимации перемещений;
]у(е) = ]у(е)(Х;У;2) = [N1 (х, у, г).. .N4 (х, у, г)] -
(3х12)-матрица базисных функций (формы
элемента); и; — (3х1)-вектор перемещений
г-го узла; и^ =[и.|...и4]т— (12Ч1)-вектор перемещений элемента.
В результате получаем уравнения для каждого элемента, которые образуют систему уравнений. Так как элементы имеют общие узлы, порядок системы может быть понижен
суммированием блоков матриц жесткости и узловых нагрузок, связанных с общими узлами. С помощью такого процесса, называемого ассемблированием, множество систем уравнений для отдельных элементов может быть объединено в одну глобальную систему уравнений для всей конечно-элементной модели [9]:
[К]{и}={Р}+{Р}+[Ь]£},
(5)
где [К] — (3пхп) — глобальная матрица жесткости конечно-элементной модели; [Ь] — (пхп)
— глобальная матрица термоупругой жесткости;
{и}= {и|...и£} — (3пх1) — глобальный вектор
узловых перемещений; {1}т ={1|..д|} — (3пх1)
— глобальный вектор узловых температур;
{р}т = {Е|Т...Р*} — (3пх1) — вектор приведенных
к узлам массовых сил; {Р}т = {Р]Т ...Р^} — (3пх 1)
— вектор приведенных к узлам поверхностных сил.
Решение системы уравнений (5) относительно перемещений имеет вид
{и}=([К]Т[К])"1 ({р}+ {Р}+М£}) (6)
Для определения радиального зазора из глобального вектора перемещений для диска и рабочей лопатки необходимо выбрать перемещения точек, находящихся на периферии этих элементов, а для статора - точек на внутренней поверхности, образующей проточную часть двигателя. Эти точки будем называть контрольными (КТ). Тогда для элемента конструкции, содержащего КТ, можно ввести матрицу связей
, состоящую из нулей и единиц, преобразующую глобальный вектор перемещений в вектор перемещений этого элемента:
А кт
{и}.
_^
;
(7)
— от действия поверхностных (газодинамических) сил:
4е-> = \у(ект) ([К]Т[К])"1 {Р}; (8)
температурные перемещения:
перемещений только для контрольных точек,
индекс (еет) можно опустить.
В работе [10] выполнен обзор существующих методов моделирования температурных перемещений и предложен способ моделирования их динамики. При моделировании зазоров основное внимание уделялось моделированию температурных деформаций, а деформации от действия силовых факторов не рассмотрены. Однако, если при помощи твердотельных моделей проанализировать общее удлинение роторных деталей турбины на максимальном режиме работы, можно сделать вывод, что вклад температурных деформаций составляет порядка 70%, а силовых факторов - 30%. Таким образом, силовые факторы вносят значительный вклад, и пренебречь их влиянием нельзя.
2. Моделирование перемещений от механических факторов
Авторы работ [11, 12] предлагают определять перемещения от действия силовых факторов следующим образом:
— от действия массовых сил:
иР =-2-иРб : пб
(10)
— от действия поверхностных (газодинамических) сил:
Р.
иР
Р Рб
г КВД6
(11)
Выражение (6) демонстрирует аддитивность перемещений, создаваемых различными факторами нагружения. Выделим следующие перемещения:
— от действия массовых сил:
и^)^е->([К]т[К])~1[Ь]£} . (9)
Так как в дальнейшем будем вести анализ
где иРб , ир^ — перемещения от действия силовых факторов, полученные на базовом режиме работы двигателя . Базовым будем считать режим работы двигателя, относительно которого определяются граничные условия теплового и силового нагружения КЭ на других режимах . Пример схемы нагружения диска приведен на рис. 1.
Из выражения (6) следует, что все перемещения, вызванные действием механического нагружения, зависят от упругих свойств материалов, которые, в свою очередь, зависят от температуры. Так как изменение температуры в диапазоне рабочих режимов значительно, ее влияние на перемещения может быть весьма существенным и, следовательно, должно учитываться при разработке моделей радиальных размеров КЭ.
Для учета этого влияния введем понятие
эквивалентной температуры 1;е , которая обеспечивает выполнение следующего условия:
|ср(х,у)-Е(х,у,1е)сЗР= |Ф(х,у)-Е(х,У)1(х,у))с1Р , Р Б
где ф(х, у) — функция координат (различная для
разных компонент тензора и факторов нагруже-
ния); Е(х,у,1;(х,у)) — модуль упругости элемента
детали с координатами (х,у) при температуре
этого элемента ^х, у); с1Р — элемент площади меридионального сечения.
Таким образом, Ц — это температура равномерно нагретого КЭ, при которой перемещение от действия силовых факторов равно перемещению при текущем распределении температурного поля КЭ и действии силовых факторов, соответствующем базовому режиму.
Рис. 1. Параметры теплового и силового нагружения диска (коэффициент теплоотдачи / температура среды)
Тогда можно записать следующую формулу подобия перемещений с учетом температурной зависимости модуля упругости:
ир
_ Е(1)-Р
иИ^) ЕОе)-рб
(12)
механического нагружения, ир(1е) — перемещение КЭ в контрольной точке только от действия механических факторов при равномерном нагреве до температуры 1;е и механическом нагружении базового режима; Е(1), — значения модуля упругости для текущей температуры в контрольной точке 1 и для эквивалентной температуры 1е ; Р, Рб — значения параметра механического нагружения на текущем и на базовом режимах. Согласно формуле
(10), массовые силы пропорциональны п2 , тогда формула для определения перемещений при текущем поле температур приймет следующий вид:
ир =
т *
(13)
где ир — перемещение КЭ в контрольной точке только от действия механических факторов при текущем поле теплового и
Е(1е) п\
Из полученной формулы следует, что для определения перемещения необходимо знать в каждый момент времени значения 1и 1:е. Данные температуры зависят от распределения температуры среды по поверхности КЭ (значения которой входят в вектор теплового нагружения) и определяются в процессе проектирования двигателя путем проведения специальных газодинамических расчетов [13, 14]. Их результаты подтверждаются или корректируются по результатам измерений температуры на участках, доступных контролю в процессе испытаний двигателя. Обычно такое распределение определяется для одного режима (базового). Пример распределения параметров теплового нагружения представлен на рис. 1.
2.1. Модель температуры среды
На произвольных режимах работы двигателя распределение не известно, но оно может быть определено по контролируемым параметрам двигателя на основе аналогии (подобия) распределения температуры торможения на произвольном и базовом режимах [15].
При выполнении обычных условий — постоянстве геометрических размеров, сохранении относительных отборов воздуха на самолетные нужды и охлаждение турбины, отсутствии влияния числа Рейнольдса на течение газовых потоков и др. — работа газотурбинного двигателя сопровождается незначительным изменением степени понижения давления турбины низкого давления и еще меньшим —
турбины высокого давления [16-18]. Изменение
* *
приведенных температур типа , на этих
режимах обычно не превышает 3 % (на рис. 2
представлены результаты расчета приведенной температуры, полученные с использованием нелинейной математической модели двух-вального ТРДД). Влияние на эти приведенные температуры условий на входе в двигатель — на порядок меньше. Это позволяет использовать распределение приведенной температуры среды на базовом режиме в качестве основы математической модели распределения температуры на других режимах.
Для диска турбины модель приведения температуры среды имеет вид
тГ
Т82 -т5*1
> 1с1ет,
(14)
где 5-1 и s2 — сечения проточной части двигателя, в которых известны значения температуры.
В предположении о постоянстве приведенных температур связь локальной температуры среды с ее значением на базовом режиме имеет вид
т*=кг(т*б-т;1б)+т;1, (15)
где индексом «б» отмечены значения величин на базовом режиме, к — коэффициент подобия режима по температуре среды, который выражается через контролируемые параметры двигателя:
к,=-
(16)
1 й2б
1в1б
Таким образом, температура среды на г-й т*/
поверхности зависит от температуры
/ Тгад
на базовом режиме и от коэффициента подобия по температуре к1.
Рис. 2. Приведенные к температуре за ТНД температуры: ^ — газа перед ТВД; ^^^ — воздуха за КВД; '4 — воздуха за КНД
2.2. Модель коэффициентов теплоотдачи
Теплопередачу на участке поверхности КЭ принято описывать критериальными уравнениями типа [14, 19]
Ыи = С-Кек-Ргт,
(17)
а-г „ р^-г „ ср'Ц .
где , , ; а, , ц,
X ц X р
— коэффициенты теплоотдачи, теплопроводности, динамической вязкости и теплоемкости при постоянном давлении; р — плотность потока; г — характерный размер; — относительная скорость; С, к, т — константы.
С помощью известных зависимостей установим связь величин, входящих в зависимость (17), с термогазодинамическими параметрами потока в проточной части двигателя.
Из уравнения состояния газа и условий газодинамического подобия следует:
р~р/Т; лу~п; Т-Т^ц^ р-р^д.
Зависимости теплофизических характеристик от температуры имеют вид [20]
0.76
! согЫ.
Подставив эти зависимости в (17), определим отношение локальных значений коэффициентов теплоотдачи на произвольном и базовом режимах — коэффициент подобия режимов по теплоотдаче:
( * > 0,8 ( * N
а п Рквд 1квд
«б * Рквд б, ,тквд б,
-0.567
(18)
Это позволяет для основных эксплуатационных режимов рекомендовать пропорциональную модель теплоотдачи:
а = аб-ка,
(19)
в которой значение коэффициента в различных условиях работы двигателя определяется формулой (18).
Таким образом, для определения граничных условий теплообмена на произвольном режиме работы двигателя необходимо задать значения локальных коэффициентов теплоотдачи на
базовом режиме аб и коэффициент подобия
режимов по теплоотдаче ка. Распределение
коэффициентов аб по поверхности детали на базовом режиме является неизменной характеристикой двигателя, определяемой в результате специальных расчетов и уточняемой по результатам испытаний двигателя. Коэффициент ка зависит как от режима работы двигателя, так и от внешних условий. Его связь с контролируемыми параметрами двигателя определяется выражением типа (18).
2.3 Модель температурного состояния КЭ
При моделировании температурного состояния с помощью МКЭ уравнение теплопроводности
dt
c-pdT"
3 к
(к
'Эх
3^ 'Эу.
д\ X
9t
3z
Зх
Зу
3z
(20)
в котором температура t = t(x,y,z,t) — неизвестная функция координат и времени,
коэффициент теплоемкости с = 1(х,у,г,1),
теплопроводности = и плотность
материала р = 1:(х,у,г,1:), можно аппроксимировать векторно-матричным дифференциальным уравнением
, (21)
ёт
где ? = — вектор неизвестных узло-
вых температур (температурного состояния
в п узлах); Т = [Т1---Т2]т — вектор теплового нагружения (температура среды); С, Л , А — матрицы теплоемкости, теплопроводности и теплоотдачи, вычисляемые на основании геометрии конечно-элементной сетки, свойств используемых материалов и граничных условий.
Для установившегося (стационарного) режима работы двигателя решение системы (21) имеет вид
tCT =(Л + А)"'-А-Т
(22)
(23)
модель №3:© зависит от уровня теплоотдачи и от коэффициента теплопроводности материала детали, изменяющегося с температурой: — к /
(эта зависимость также аппроксимировалась полиномом).
В качестве исходных данных для формирования этих моделей использованы значения температуры , полученные с использованием моделей высокого уровня. В качестве внешнего теплового воздействия задавались значения параметров Т8*1, кг и ка . Значения 0i для каждого 1-го сочетания этих параметров были рассчитаны на основании формулы (23). Значения ©Р были получены для тех же условий с помощью модели высокого уровня и принимались в качестве эталонных. Тогда относительное отклонение характеризует погреш-©? _
ность модели параметра 0, а абсолютное отклонение Д! = ^, - ^} характеризует погрешность модели температуры. Результаты расчета этих погрешностей приведены в табл. 1.
Таблица 1
Погрешности моделей определения 0 и ^
Формула (15) преобразует алгоритм расчета температуры в контрольной точке на установившемся режиме из операций с конечно-элементными матричными моделями (21) и (22) в простейшую алгебраическую операцию с зависимостью приведенной температуры от режима работы и с двумя контролируемыми параметрами двигателя:
8©, % At, К
модель №1 ~2% ~5
модель №2 ~0,5 ~1,5
модель №3 ~0,04 ~0,2
где © — коэффициент пропорциональности, который можно рассматривать как приведенное к контролируемым параметрам значение температуры в критической точке. Этот коэффициент зависит от режима работы двигателя.
Для формирования этой зависимости рассмотрено три альтернативных модели:
модель №1: 0 = const (эта константа была определена для базового режима);
модель №2: 0 зависит от уровня теплоотдачи:
0 = f(ka) (данная зависимость аппроксимировалась полиномом);
Анализ приведенных выше значений погрешностей дает основания в дальнейшем использовать модель ® = которая
обладает наибольшей точностью. В этой модели к^ =—— параметр, характеризующий ЦТб)
отличие коэффициента теплопроводности материала от его значения на базовом режиме;
Т, Тб — характерные (для учета изменения теплопроводности материала) температуры на рассматриваемом и базовом режимах. Характерной для всего узла конструкции считалась температура, равная средней по поверхности температуре среды, взвешенной по локальному коэффициенту теплоотдачи: |аТ*ёР
JadF
где ей7 — элемент поверхности.
Характерную температуру Т можно найти, обработав граничные условия базового режима с использованием формулы (24). Далее можно найти ее приведенное значение
. (25)
Т52б -Т81б
В дальнейшем для всех установившихся режимов характерную температуру можно рассчитывать, не используя высокоуровневые модели, по следующей формуле:
т=0-(т;2-т;1)+т;1. (26)
— к /
Для формирования зависимости при помощи твердотельной конечно-элементной модели (модели верхнего уровня) для установившихся режимов работы двигателя был рассчитан ряд температур ^. Условия, при которых были рассчитаны эти температуры, отличались значениями Т8* , коэффициентами подобия по температуре к<. (см. формулу (15)) и теплоотдачи ка , которые покрывают все возможные сочетания этих коэффициентов, которые могут встретиться при работе двигателя (даже самые редкие сочетания). Из формулы (23) были определены значения ©. Для этих же сочетаний параметров теплового нагружения по формуле (26) с использованием ранее полученного коэффициента© для базового режима по формуле (25) была определена характерная температура Т;. С использованием зависимости коэффициента теплопроводности материала от температуры были получены Х(Тб) и Х(Т), а затем рассчитан коэффициент к^ на этих режимах. Полученная зависимость монотонна, однозначна и хорошо аппроксимируется полиномом четвертой — пятой степени (рис. 3).
Сравнение значений ^ , полученных с использованием моделей высокого уровня, со значениями, полученными по упрощенной модели с использованием формулы (23), дает максимальную абсолютную погрешность 0,32 К, что соответствует относительной погрешности около 0,05%.
О 0018
О -
О 0.5 1 1.5 2
Рис. 3. Результаты расчета © : ® — смоделированные значения: --аппроксимирующая функция
Температурное состояние КЭ на неустановившихся режимах работы двигателя можно определить, решая систему уравнений (21) при начальном состоянии в контрольной точке
t(T0) = t0. Решение можно представить в следующем виде:
т
tco = *0 + J nt (t. л) ■ dtCT (Л), (27)
где tCT(ri) — стационарная температура контрольной точки при условиях, существующих в момент времени h, рассчитанная по формуле (23); nt(T, ti) — переходная характеристика температуры КЭ при ступенчатом изменении
tCT в момент времени ti. Эта формула позволяет заменить алгоритм расчета температуры в контрольной точке при помощи МКЭ на операции с переходными характеристиками . Основная проблема создания динамической модели (модели на переходных режимах) вида (27) — необходимость расчета переходных характеристик для всевозможных вариантов изменения граничных условий теплообмена.
В процессе разработки моделей на неустановившихся режимах были получены переходные характеристики при ступенчатом изменении коэффициентов теплоотдачи. Расчеты проводились с использованием моделей верхнего уровня.
В результате были получены кривые переходных процессов температуры, для которых использовалась следующая система обозначений:
t(T) = t(T, к« к£>), где к® — начальное значение параметра теплоотдачи; к£Р — измененное (текущее) значение параметра теплоотдачи.
По результатам определялись нормированные переходные характеристики температуры при ступенчатом изменении коэффициента теплоотдачи:
, (28)
*ст_1:0
где ^ и 1ст — начальное и конечное (при т —> оо) значения температуры.
Расчет выполнялся при изменении уровня ка в момент времени т = 0 с к® =0 на величину . Температура среды по поверхности диска также изменялась скачком. Начальное распределение температуры по всем участкам было равно 293 К, а в момент времени т = 0 изменялось до базового распределения. Результаты представлены на рис. 4.
Исходя из характера переходных тепловых процессов, было предложено описывать их апериодическими звеньями высокого порядка с переменными коэффициентами, передаточные функции которых имеют вид
Х(8) + 1 '
(29)
где — выходной параметр, в данном случае температура КЭ 1(т)в момент времени т,
Х(в) — воздействие, в данном случае 1ст, т — количество экспонент, при помощи которых описывается переходная характеристика,
к^ка), Т^ка) — коэффициент усиления и
постоянная времени ] -ой экспоненты для текущего режима.
1 1 , т. ...
V 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
50
100
150
200
250
Рис. 4. Ступенчатые переходные характеристики температуры диска:
к£} = 0,2 - к£} = 0,4 ; - к£} = 0,6 ;
^-к«=0,8;—-к£>=1; -к« =1,2.
Как видно из рис. 5, постоянная времени переходных процессов в наибольшей степени зависит от интенсивности теплоотдачи, которая, в свою очередь, наиболее точно характеризуется коэффициентом подобия по теплоотдаче для данной детали. Тепловое состояние (температура) диска с достаточной точностью описывается апериодическим звеном третьего порядка, коэффициенты которого зависят от режима работы двигателя, представленного значением ка. Эти зависимости представлены на рис. 5 и 6.
При программной реализации уравнения (27) интеграл заменяется на сумму переходных характеристик. Изменение температуры 1ст представлено как последовательность приращений Д1ст 1. В результате время расчета одной точки значительно увеличивается, что не позволяет моделировать процесс в реальном масштабе времени. Для решения этой проблемы уравнение (29) представлено в дискретной форме, тогда переходный процесс принимает вид
Уу = ' +
Ъ ГХ;
и
(30)
(31)
где у^ — значение выходного параметра, рассчитанное для j -й экспоненты (всего т экспонент) при 1 -х параметрах теплового нагружения; — значение выходного параметра, рассчитанное для ] -й экспоненты при параметрах теплового нагружения, соответствующихпредыдущему шагу 1
интегрирования; - —; bj = 1 - aj — констан-
— + 1
Т-
ты .¡-й экспоненты; - шаг дискретизации; Ту , ку — постоянная времени и коэффициент
усиления ] -й экспоненты при 1 -х параметрах теплового нагружения.
т. с
а
Рис. 7. Погрешность температурной модели на неустановившихся режимах
2.4 Определение эквивалентной температуры
Применив формулу (12) к условиям базового режима, получим выражение для определения значения эквивалентной температуры на базовом режиме:
Рис. 5. Постоянные времени экспонент, при помощи которых описываются переходные процессы, изображенные на рис. 4: а — первая экспонента; б — вторая экспонента; в — третья экспонента
На рис. 7 приведены результаты расчета температуры при помощи моделей высокого уровня и при помощи выражения (31). Данный переходный процесс имитировал запуск двигателя при температуре окружающей среды +20 0С, выход на режим малого газа с выдержкой на нем в течение 120 с, затем выход на базовый режим. При этом параметры теплового нагру-жения изменялись скачкообразно.
Максимальное отличие результатов наблюдается в момент изменения режима и составляет 3-7 К. Эту погрешность, существующую кратковременно, можно существенно уменьшить за счет уменьшения шага дискретизации АХ. На рис. 7 приведены результаты для Д1=0,5 с. При уменьшении А! до 0,1 с максимальное расхождение уменьшилось до 4 К.
up (tee)
Рис. 6. Коэффициенты усиления экспонент, при помощи которых описываются переходные процессы, изображенные на рис.4: — первая экспонента; — вторая экспонента;
^^^ — третья экспонента
(32)
иРб Е(*б) Для определения ^ предварительно с использованием модели высокого уровня были рассчитаны перемещения диска от действия центробежных сил ир(1е), значения которых соответствуют базовому режиму, при равномерном нагреве до различных температур 1£. Диапазон варьирования температуры должен покрывать рабочий диапазон двигателя. Далее по модели высокого уровня необходимо определить перемещения иРб при базовом распределении температур и базовом механическом нагружении. Температура в контрольной точке диска на этом же режиме определяется по формуле (23).
Далее для всех температур необходимо
рассчитать отношения перемещений
Цр(*е) .
Используя изменение механических свойств материала по температуре, определяется от-
ношение
E(te) E(t6)
модуля упругости материала
для температур te, при которых определены значения up(te), к модулю упругости при базовой температуре детали t6 . Точка пересечения этих двух графиков дает искомую температуру
teg (см. рис. 8).
в
Рис. 8. Графическое определение эквивалентной температуры на базовом режиме ^ :
Е(1е)
Е0б) ирб
Для определения эквивалентной температуры на произвольном режиме воспользуемся следующим следствием подобия температурных полей, справедливым в широком диапазоне изменения ка (условий теплообмена):
1е-Т Ц. ~тб 1-Т С-Тй
(33)
где I - значение температуры в контрольной точке для данного режима, рассчитанное по формуле для установившегося режима работы двигателя, либо по формуле (27) для переходных режимов;
1;б - значение температуры в контрольной точке для базового режима; Т , Тб - характерная температура на произвольном и на базовом режимах работы двигателя (см . формулу (26));
^ - эквивалентная температура на базовом режиме, получаемая в результате решения уравнения (32); 1:е - эквивалентная температура на произвольном режиме.
Определив один раз значение т] для условий базового режима, будем рассчитывать эквивалентную температуру по формуле
3. Верификация модели перемещений
(34)
Для проверки корректности алгоритма моделирования деформации диска от действия центробежных сил с использованием модели высокого уровня были определены значения деформации для различных сочетаний параметров теплового нагружения , ка ,
и
этих же сочетаний параметров по формуле (13)
определены перемещения иР только от силовых факторов с учетом влияния температуры на модуль упругости.
Результаты сравнения погрешности определения радиальных перемещений, полученных с использованием формулы (13) (т.е. с учетом зависимости модуля упругости от температуры), и модели (10) (т.е. без учета этой зависимости), представлены на рис. 9. Анализ показал, что погрешность модели с учетом зависимости модуля упругости от температуры не превышает 1,61 %, в то время как погрешность модели без учета изменения модуля упругости достигает 20 %. Следовательно, модель (13) можно рекомендовать к использованию при моделировании радиальных размеров КЭ, т.к. она дает удовлетворительную погрешность.
На рис. 10 приведены результата расчета суммарных перемещений диска от действия силовых факторов, а также температурные деформации, рассчитанные с использованием модели высокого уровня и с использованием формулы (6). В качестве исходных данных задавалось совместное изменение теплового и силового нагружений. Характер изменения внешних условий следующий:
В период времени от 0 до 120 с - ка = 0.208,
к4 = 0.4857 , т;, = т£вДмг , кп = 0.6654 ;
В период времени от 120 с до 500 с - ка = 1,
]Ч=1, ^ == ^КВДб , кп=1.
Максимальное отличие результатов наблюдается в момент изменения режимов и составляет 0,11 мм, далее это расхождение уменьшается.
/ /
1 ■ I ■ ■ Г
■ ■ ■ ■
1.1% ■а- ■ ч—— -V ■ ■
-1.61% А*— — ■ -«Н* л!-
параметра силового нагружения кп =
. Для
Л,
Рис. 9. Погрешности моделей определения перемещений диска от действия центробежных сил на установившихся режимах: ♦ - модель с учетом изменения модуля упругости от
температуры; - модель без учета модуля упругости от температуры
пб
1 J__
| 0.07mv to 02мм
У
'to 11 м» i
О 100 200 300 400 L- L
Рис. 10. Погрешность определения суммарных перемещений диска
Заключение
В результате анализа температур и деформаций диска получены основные зависимости между его радиальными деформациями и параметрами проточной части двигателя. Предложен простой метод учета температурного поля детали при расчете перемещений от действия силовых факторов, т.е. учета изменения модуля упругости от температуры.
Полученные результаты позволили рекомендовать модель (31) для определения температуры диска на неустановившихся режимах и модель (13) для расчета перемещений от действия силовых факторов.
Данный подход позволяет анализировать динамику перемещений для любой детали двигателя, подверженной механическому нагруже-нию и воздействию высоких температур.
Литература
1. Sobey A.J. Control of aircraft and missile powerplants [Text] / A.J. Sobey, A.M Suggs. - New York: Wiley, 1963. - 432p.
2. Гуревич О.С. Интегрированное управление силовой установкой многорежимного самолета [Текст] / О.С. Гуревич, Ф.Д. Гольберг, О.Д. Селиванов; под общей ред. О.С. Гуревича. — М.: Машиностроение, 1993. - 304с.
3. Saravanamuttoo H.I.H. Simulation of gas turbine dynamic performance [Электронный ресурс] / H.I.H. Saravanamuttoo, A.J. Fawke // Proceedings of the ASME. Gas turbine conference, Brussels, Belgium, 1970. — ASME Paper No 70-GT-23 — 1 электрон.опт.диск (CD-ROM).
4. Jaw L.C. Aircraft engines controls: design, system analysis, and health monitoring [Электронный ресурс] / L.C. Jaw, J.D. Mattingly // American Institute of Aeronautics and Astronautics, Inc, Reston, Virginia, 2009. — 1 электрон.опт.диск (CD-ROM).
5. Kurzke J. Transient simulations during preliminary conceptual engine design [Электронный ресурс] / J. Kurzke. — 2011, ISABE-2011-1321. — 1 электрон.опт.диск (CD-ROM).
6. Добрянский В.Г. Динамика авиационных ГТД [Текст] / В.Г. Добрянский, Т.С. Мартьянова. — М.: Машиностроение, 1989. — 240 c.
7. GasTurb 12 A program to calculate design and off-design performance of gas turbines. Users manual [Text] / GasTurb GmbH: Dr.-Ing. Joachim Kurzke. — Aachen, 2012. — 329 p.
8. Олейник А.В. Эксплуатационный мониторинг температурного состояния детали газотурбинного двигателя как задача динамики конечно-элементной модели в пространстве состояний [Текст] / А.В. Олейник / Авиационно-космическая техника и технология. - 2004. - №4/12. — С. 38-42.
9. Еременко С.Ю. Методы конечных элементов в механике деформируемых твердых тел [Текст] / С.Ю. Еременко. — X.: Основа, 1991. — 272 с.
10. Yepifanov S.V. Modeling the GTE under its dynamic heating conditions [Text] / S.V. Yepifanov, R.L. Zelenskyi, I.I. Loboda // Engineering for gas turbines and power (March 2015). -2015. Vol.137. —Number 3.-P.031506-1-031506-10.
11. Kypuros J.A. A reduced model for prediction of thermal and rotational effects on turbine tip clearance [Text] / J.A. Kypuros, K.J. Melcher. — Tech. Rep. NASA, TM-2003-212226, 2003 (http://ntrs.nasa.gov/archive/nasa/ casi.ntrs.nasa.gov/20030032933.pdf).
12. The whole-engine model for clearance evaluation [Электронный ресурс] / A.N. Arkh-ipov, V.V. Karaban, I.V. Putchkov, G. Filkorn [at al.] // Proceedings of ASME Turbo Expo 2009, 8 —12 June 2009.— Orlando, Florida, USA, 2009. — ASME Paper GT2009-59259. — 1 электрон.опт.диск (CD-ROM).
13. Копелев С.З. Конструкции и расчет систем охлаждения ГТД [Текст] / С.З. Копелев, А.Ф. Слитенко. — X.: Основа, 1994. — 240 с.
14. Швец И.Т. Воздушное охлаждение деталей газовых турбин [Текст] / И.Т. Швец, Е.П. Дыбан. — К.: Наук. думка, 1974. — 487 с.
15. Олейник А.В. Концепция и методы мониторинга выработки ресурса авиационных ГТД на основе идентификации динамики температурного и напряженного состояния основных деталей: дис д-ра техн. наук: 05.07.05; защита 17.10.06; утв. 18.03.07г. / Олейник Алексей Васильевич. — X, 2006.- 240 с.
16. Кулагин В.В. Теория газотурбинных двигателей [Текст] / В.В. Кулагин. Кн. 2. Совместная работа узлов, характеристики и газодинамическая доводка выполненного ГТД.— М.: МАИ, 1994. — 304 с.
17. Шляхтенко С.М. Теория воздушно-реактивных двигателей [Текст] / С.М. Шляхтенко. — М.: Машиностроение, 1975. — 568 с.
18. Шляхтенко С.М. Теория двухконтурных реактивных двигателей [Текст] / С.М. Шляхтенко, В.А. Сосунова. - М.: Машиностроение, 1979. - 432 с.
19. Кутателадзе С.С. Анализ подобия в те-плолфизике [Текст] / С.С. Кутателадзе. - Новосибирск: Наука, 1982. - 280 с.
20. Копелев С.З. Тепловое состояние элементов конструкции авиационных двигателей [Текст] / С.З. Копелев, С.В. Гуров. - М.: Машиностроение, 1978. - 208 с.
Поступила в редакцию 12.05.2015
Р.Л.Зеленський, С.В. Сшфанов, С.В. Марценюк, В.В. Бойко. Моделювання динамши рад1ального розм1ру диска турбши, зумовлено! д1ею мехашчних фактор1в
Важливим параметром, який впливае на ККД турбти, е рад1альний зазор М1Ж ротором i статором. Розвиток систем активного керування зазором зумовлюе необх^дшсть розробки i використання математичних моделей, ят вiдбивають динамку змти зазору в рiзних умовах роботи двигуна. Значения радiального зазору визначаеться радiальними розмiрами деталей, я^ приймають участь у його формуванш: диска, лопатки i корпуса. У цш роботi розглянуто моделювання радiального розмiру диска. Для цього необхiдно ураховувати не тыьки нагрiв, але i дт силових факторiв навантаження. Використання спрощених моделей вiдповiдних деформацт не дае прийнятног точностi, тому що спро-щена модель не ураховуе змти мехашчних властивостей (зокрема модуля пружност1 матерiалу) вiд температури. Запропоновано методику визначеннярадiального подовження диска вiд дИ силових факторiв iз динамiчною корекщею значення модуля пружностi за температурою. Порiвняльнiрезультати моделювання подовження iз використанням МКЕ i спрощеног моделi з корекщею модуля пружностi показали, що запропонована модель забезпечуе необхiдну точшсть.
Ключов1 слова: динамiчна математична модель, радiальне подовження, мехашчне навантаження, температура, перехiднi характеристики, граничш умови.
R.L. Zelenskyi, S.V. Yepifanov, Ye.V.Martceniuk, V.V.Boiko. Simulation of a turbine disk radial size under the mechanical loading
The significant parameter that acts turbine efficiency is the radial clearance between rotor and stator. Development of the radial clearance active control systems makes necessity of mathematical models development and application, which represent dynamics of clearance variation in different engine operational conditions. The value of radial clearance is determined by radial size of parts that form it: disc, blade and casing. This paper represents simulation of the disc radial size. Action of power factors must be added to action of heating. Simplified models of the respective strains not provide required precision, as the simplified model not takes into account variation of mechanical properties (particularly modulus of elasticity) with temperature. Method of the disc elongation caused by power factors is proposed. This method includes correction of elasticity modulus with temperature. Comparative simulation of elongation using finite-element and simplified models showed that the proposed model provides required precision.
Key words: dynamic mathematical model, radial elongation, mechanical loading, temperature, transient performances, boundary conditions.
