КОНТАКТНАЯ ЗАДАЧА В АНАЛИЗЕ ДИНАМИЧЕСКОГО ПОВЕДЕНИЯ СБОРНЫХ РОТОРОВ ТУРБОМАШИН Текст научной статьи по специальности «Механика и машиностроение»

ISSN 1814-1196

http://journals. nstu. ru/vestnik Science Bulletin of the NSTU Vol. 60, No. 3, 2015, pp. 113-129

Научный вестник НГТУ том 60, № 3, 2015, с. 113-129

моделирование процессов и устройств

modeling of processes and devices

УДК 539.3/.6

Контактная задача в анализе динамического поведения сборных роторов турбомашин

М.А. ДУДАЕВ1, А.А. ПЫХАЛОВ2

1 664074, РФ, г. Иркутск, ул. Лермонтова, 83, Иркутский национальный исследовательский технический университет, аспирант. Е-таИ: [email protected]

2 664074, РФ, г. Иркутск, ул. Лермонтова, 83, Иркутский национальный исследовательский технический университет, доктор технических наук, профессор. Е-тай: [email protected]

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

Ключевые слова: сборный ротор, контактная задача, метод конечных элементов, математическая модель, критическая скорость вращения, динамическое поведение, метод Ньюмар-ка, переходный процесс, колебания

DOI: 10.17212/1814-1196-2015-3-113-129

Статья получена 30 апреля 2015 г.

ВВЕДЕНИЕ

Современное развитие роторных систем турбомашин основывается на принципах многокаскадности и повышения рабочих параметров [1-3]. При этом основополагающим подходом в определении динамических параметров ротора, в частности критических скоростей (частот) его вращения, остается натурный эксперимент, или гонка его реальных прототипов. Одной из главных причин такого подхода является то, что роторы большинства высоко-нагруженных турбомашин являются принципиально сборной конструкцией, жесткость которой из-за наличия контактных сопряжений в значительной степени отличается от жесткости монолитного аналога и изменяется в ходе раскрутки ротора. Таким образом, неизвестным остается интерпретация влияния условий сопряжения деталей на динамические параметры роторов современных высоконагруженных и многокаскадных турбомашин. Например, это относится к авиационным газотурбинным двигателям [4, 5] как к наиболее энергоемким механическим системам. Сборная конструкция их роторов отличается большим количеством сопрягаемых деталей и разнообразием конструктивных решений по их сопряжению.

Для решения представленной проблемы требуется проведение подробного инженерного анализа, в особенности изменения критических частот вращения ротора [6], а также учет ряда физических процессов механики деформируемого твердого тела. В частности, к таким процессам относится контактное взаимодействие деталей ротора, гироскопический эффект, влияние на ротор поля центробежных сил и другие. Эти эффекты способствуют изменению критической частоты вращения реального сборного ротора по сравнению с математической моделью его монолитного аналога. Наиболее актуальным в этом случае является изучение опасного приближения критической частоты вращения к диапазону рабочих частот ротора (рис. 1), а также повышение общего уровня его амплитуд колебаний.

0.16 0.14

I 0.12

1 0.10

ч>

| 0.08 8" 0-06 ^ 0.04 0.02

0 50 100 150 200 250 300 350 400 450 п, об/с

Рис. 1. Амплитудно-частотная характеристика роторов

Рабочий диапазон Передняя опора

Модель сборной конструкции ротора Модель ротора с упругими элементами Модель монолитного аналога сборного ротора

Определенной альтернативой натурному эксперименту (в отношении сокращения материальных и временных затрат при проектировании и доводке конструкции сборного ротора, а также для повышения уровня информативности о его работоспособности) служит развитие теоретических подходов и, соответственно, расчетных методов анализа динамики сборных роторных систем. На сегодняшний день наиболее эффективным из них является метод конечных элементов (МКЭ), использование которого построено с реализацией решения контактной задачи теории упругости [7, 8], позволяющей моделировать и, соответственно, изучать особенности работы тех или иных конструктивных решений сборного ротора, т. е. контактное взаимодействие или условия сопряжения его деталей. Ниже представлена математическая модель и определенные результаты предлагаемого для решения показанной выше проблемы.

1. МАТЕМАТИЧЕСКАЯ МОДЕЛЬ

Основное уравнение движения упругой системы под действием внешних сил с применением МКЭ имеет вид [7, 8]

[М]{5} + [С]{5} + [^]{5} = {^}, (1)

где [М ] - матрица масс, изменяющаяся с течением времени ввиду изменений инерционных характеристик ротора в ходе раскрутки; [С ] - матрица коэффициентов эквивалентного вязкого демпфирования; [А"] - матрица жесткости сборной конструкции ротора; , |<51 и ¡8} - соответственно векторы ускорений, скоростей и перемещений узлов конечно-элементной модели, являющиеся основными неизвестными; - вектор внешних сил, приложенных к узлам конечно-элементной модели и являющихся функциями времени.

Матрица жесткости сборного ротора [К ] по сравнению с матрицей

жесткости монолитного аналога имеет более сложную структуру. Главными ее составляющими являются матрица конструкционной жесткости деталей ротора [ К ] и матрица контактной жесткости [ К ] , используемая для математического моделирования сопряжений деталей. Кроме того, в ней учитывается матрица жесткости от гироскопического эффекта ротора [К] и матрица жесткости от начальных напряжений ротора [ К] :

[КНК]к +[К]конт +[К]г +[КI. (2)

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

Вектор узловых сил включает в себя силы от статического и динамического дисбаланса ротора {-^}д, силы контактного взаимодействия деталей

ротора {^}конт, центробежные {-^}ц и гравитационные силы :

И = Ид + П«ЖГ +{^}ц + Пр • (3)

Рассмотрим более подробно компоненты матрицы жесткости согласно уравнению (2).

Матрица конструкционной жесткости [ К ] представляет собой структуру из симметричных, положительно определенных матриц деталей ротора, каждая из которых формируется суммированием из матриц конечных элементов (КЭ) и определяется выражениями [9, 10]

[ К ]к =

[ к I

[ к ]

[ к ],

(4а)

Ые

[К] = ![КГ. ;

у=1 г.

[к]е, = 11 щ [,

V

(4б)

(4в)

где [ К ] - матрица конструкционной жесткости 1-й детали; Ые - количество КЭ, составляющих 7-ю деталь; [К]е - матрица жесткости у-го КЭ в детали;

|В| - матрица градиентов функций форм КЭ в полярно-цилиндрической системе координат, переводящая вектор узловых перемещений КЭ и = {иг ид иг} в вектор деформаций в = {вг Вд вz угд УQz угг} [9, 10]; [Б] - матрица упругих постоянных материала КЭ [11, 12]; [V] - объем КЭ.

Отдельного рассмотрения требует вопрос формирования матрицы контактной жесткости [К]конт, совместно с которой формируется вектор контактных сил {^}конт. Сопряжения деталей ротора моделируются специальным типом конечного элемента (рис. 2), который в работах [7, 8] назван контактным элементом сопряжения конструкций - КЭСК; при этом для моделирования сопрягаемых деталей используется подход узел-в-узел.

Рис. 2. Моделирование контактного взаимодействия деталей: с зазором (слева) и с натягом (справа)

Исходное состояние контактного элемента характеризуется величиной его начального раскрытия , устанавливаемого на этапе моделирования.

При моделировании натяга Мо < 0 начальное состояние КЭСК - «закрыт», при моделировании зазора Мо > 0 состояние КЭСК - «открыт». В случае моделирования сопряжения по номинальному размеру Мо = 0 КЭСК - «закрыт». В ходе решения системы уравнений МКЭ состояние каждого КЭСК может изменяться от состояния «открыт» до состояния «закрыт» и наоборот, в связи с чем решение задачи проводится итерационно и останавливается на к-й итерации в случае, если состояние большинства или всех КЭСК не изменилось по сравнению с итерацией к - 1.

Оценка состояния КЭСК производится на основании анализа невязки (неравенства) поля перемещений его узлов:

АМ = МЛ - МБ, + % , (5)

где и л и и в - перемещения узлов А и Д-го КЭСК (рис. 2).

При этом КЭСК «закрыт», если Аиг- < 0, и «открыт», если Аиг- > 0. Смена состояний КЭСК в ходе решения обусловливает модификацию глобальной матрицы жесткости [ К ] прямо в ходе решения за счет добавления в нее матрицы контактной жесткости:

^КЭСК

[К] = У <

У -1конт

г=1

[К ]кэск , ам > 0

[ К ]КЭСК' АМ- <

где матрица жесткости открытого КЭСК

[ * ]

+

КЭСК

[*]+ -[*]+ -[*]+ [*]+

, [*]+=

* „

к

(7)

- закрытого КЭСК: [ * ]КЭСК

[к]- -[*]--[*]- [*]-

, [*]-=

* -

к-

*-

(8)

где К+ и К+ - жесткость открытого КЭСК в осевом и тангенциальном направлениях соответственно - малая величина, необходимая и достаточная для обусловленности глобальной матрицы жесткости - величина, назначаемая согласно [7, 8] на два-три порядка меньше минимальной жесткости в узлах КЭ детали; Кп и Кт - жесткость (штрафная) закрытого КЭСК соответственно в осевом и тангенциальном направлениях - величина, подбираемая согласно [7, 8] на два-три порядка больше максимальной жесткости в узлах КЭ детали.

Вектор контактных сил вычисляется в соответствии с выражениями

{F}

^ЮСК |о, Ащ > 0, i=1 1 {F}контг ,

AUj < 0,

(9а)

{^конт/ = *~-'Аи I1 00 -1 0 °f. (9б)

Контактные силы, приложенные к узлам КЭСК, имеют противоположные знаки и устраняют натяг закрытого КЭСК.

Процедура решения контактной задачи показана на рис. 3 на примере контакта двух пластинчатых элементов. На начальном этапе решения задачи пластины располагаются одна над другой с зазором U0 > 0 ; перемещения узлов ид = Uß = 0. По формулам (9) имеем Au > 0 , следовательно, контактный

элемент открыт. Жесткость контактного элемента [* ] = r L -1конт

-3

= min(*i, *2) • 10 , контактных сил нет. После решения задачи перемещения узлов ид > Uß - U0 невязка поля перемещений КЭСК Au < 0, что говорит об изменении состояния КЭСК на закрытое. В свою очередь, само по себе изменение состояния КЭСК говорит о необходимости следующего шага решения задачи. Теперь КЭСК считается закрытым, жесткость КЭСК

[*]конг = max(*1, *2)-103, модули контактных сил {^}конг =[*]конг •Н .

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

Состояние КЭСК Состояние КЭСК

меняется в ходе расчета не меняется в ходе расчета

Рис. 3. Укрупненный алгоритм решения контактной задачи

Матрица дополнительной жесткости, обусловленной гироскопическим моментом ротора [7, 13], имеет вид

[ * !г =О2

2 1 О ,

т;

(10)

где О - угловая скорость прецессии ротора; ю - угловая скорость вращения ротора; т1 - масса, сосредоточенная в 7-м узле ротора от сопряженных в этом узле КЭ.

Матрица дополнительной жесткости от начальных напряжений [9, 14] определяется выражениями

Ые

[ * 1 = 1 [ * е ;

г=1

[ * Е=1 К н\о\ж,

(11а) (11б)

где [К- матрица жесткости от начальных напряжений /-го элемента; . - матрица производных функций форм N конечного элемента в поляр-

но-цилиндрической системе координат:

дг

( М 1 дМ > — +--1

V г г д9 дг

=

I ч

дг

( 1 дК ^

V г

г де

дК

дг

«I

(12)

п - количество узлов конечного элемента; I- единичная матрица (3^3); [а]. -матрица, составленная из компонент напряжений КЭ:

н=

С г1 Тге1 ^гг1

Се1 тег1

Бут

С г1

(13)

Рассмотрим силы, возбуждающие колебания ротора. Основными из них являются силы и моменты от дисбалансов ротора [5, 8]: статического Ое = те, где т - масса ротора, е - радиальный эксцентриситет, и динамического = I% , где I - момент инерции ротора относительно оси симметрии; X - угловой эксцентриситет. Статический дисбаланс возникает ввиду смещения центра масс ротора от оси вращения (рис. 4, а), динамический - ввиду несовпадения главной оси инерции ротора с осью вращения (рис. 4, б).

Центр масс

1 'лавная ось инерции

Ось вращения

/

Ось вращения а б

Рис. 4. Дисбалансы ротора

Для задания дисбалансов в модели используется специальный экстраполирующий конечный элемент (рис. 5).

Рис. 5. Моделирование дисбалансов

Указанный элемент соединяет узлы КЭ модели диска ротора с реперным узлом, расположенным строго на оси вращения. Для реперного узла задается статическая и динамическая неуравновешенность, которая пересчитывается на силы, приложенные к узлам диска в соответствии с выражением [8] вида

Ид =

га2т

-е, (8т ф e 008 е,- + 008 ф e 8т е,) ei (8Ш ф е 8Ш е, - 008 ф е 008 е,)

х, • г (8т ф % 008 е, + 008 ф % 8т е,)

(14)

где т1 - масса, сосредоточенная в /-м узле диска от связанных с ним КЭ; е1 и х, - эксцентриситеты /-го узла диска:

е, = еЯ +(г,-гя)Хя, X, = Хя,

(15)

ф е и ф% - угловые координаты положения эксцентриситетов ея и хя; г, и е, - полярные координаты /-го узла диска.

Решение уравнения (1) сводится к его численному интегрированию с применением метода Ньюмарка [7, 8, 15]. Основная идея метода заключается в разделении исследуемого временного интервала переходного режима ротора (Т) на малые конечные промежутки времени & (рис. 6, а). В таком случае, используя принцип Даламбера, условия динамического равновесия системы строго выполняются лишь в определенные моменты времени , что приводит уравнение (1) к виду

(16)

где п - количество временных интервалов At.

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

зовать 12-16 шагов интегрирования на каждый период. В этом случае выгоднее использовать адаптивный временной интервал Д^- (рис. 6, б).

Рис. 6. Постоянный (а) и адаптивный (б) временные шаги интегрирования

Согласно методу Ньюмарка, ускорение в пределах каждого временного отрезка усредняется и считается постоянной величиной (рис. 7).

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

[*]дг I5}- , (17)

где [ К ] - динамическая матрица жесткости:

К1д, =№—М]; (18)

- динамический вектор сил:

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

МЛ^срГ^' (20)

где [М] - приведенная матрица масс; - приведенный вектор сил.

Выражения для [M ] и } рассматриваются ниже.

Vil °ср i т 5г

t, 1 - - м- _

ti

5,-1 t 5,

Рис. 7. Аппроксимация параметров движения по методу Ньюмарка

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

щ

^ >с

H={S}m+HmAí + M

A ti

срг 2

(21а) (21 б)

(21 в)

Величины с индексом / - 1 считаются известными, поскольку вычислены на предыдущем шаге интегрирования. Если среднее ускорение текущего шага вычислено, могут быть рассчитаны и ускорение , скорость {§}

и перемещение {5}. текущего шага интегрирования.

Подстановка указанных величин в уравнение движения приводит его к

виду

+

б

б

5

срг 2

(22)

После раскрытия скобок имеем

2ИЧСР "МП,+"1СШ

Фг

Л,2

Группировка известных значений справа и неизвестных слева дает

(23)

А/2

2[ М ] + А/ [С ].+ —[ К ]

N

+ [^]{8}м -Щ + ^ М]{5}м - [*]{8}м .

(24)

Таким образом, полученное уравнение имеет вид выражения (20), в котором приведенная матрица масс определяется по формуле

[МI =

А/2

2[ М ] + А/ [С ].+—[ К ]

(25)

а приведенный вектор сил

Ч^ + МН,! "[И + А,'"[<1Н-1 - [^]{8}м- (26)

На основании изложенного подхода проведена программная реализация решателя для динамического анализа сборных роторов. Анализ сходимости решения проводился для каждого типа КЭ и физических процессов, речь о которых шла выше, в отдельности.

2. ЧИСЛЕННЫЙ ЭКСПЕРИМЕНТ

В качестве тестовой расчетной модели рассматривается конструкция сборного ротора, показанного на рис. 8.

Ротор установлен на упругие опоры с жесткостью порядка 105 Н/мм, что соответствует жесткостям реальных подшипниковых опор [7]. Диск ротора соединяется с валами посредством фланцев. Моделирование стяжных болтов осуществляется с применением КЭ балочного типа. Применение такого подхода обосновано в работе [7]. Поскольку болт в соединении подвергается совместному действию изгиба и сдвига, используется балочная модель Тимошенко [16]. Величина радиального эксцентриситета составляет 0,05 мм, а углового 0,001 рад. Масса ротора в сборе составляет 232 кг.

Болтовое соединение (балочный КЭ)

Поверхности сопряжения

Рис. 8. Тестовая модель сборного ротора:

а - эскиз ротора; б - схема закрепления; в - сопряжения ротора

Рис. 9. Конечно-элементная модель сборного ротора (а), низшие формы

колебаний ротора (б)

КЭ модель ротора представлена на рис. 9, а. Ротор раскручивается от нуля до 400 об/с с постоянным угловым ускорением 200 об/с2. Общее время раскрутки составляет 2 с, а число полных оборотов, совершенных за время раскрутки, равно 400. Подбор расчетного ускорения подробно рассмотрен в работе [8]. Снижение углового ускорения приводит к избыточному количеству

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

В ходе раскрутки ротор проходит три низшие критические скорости. Им соответствуют формы колебаний, представленные на рис. 9, б. Из них первые две формы обусловливаются податливостью опор, в связи с чем называются «твердотельными» формами колебаний, а третья обусловлена дополнительной податливостью самого вала, в связи с чем называется «изгибной» формой.

В результате расчета получены графики зависимости поступательных перемещений на упругих опорах ротора, представленные на рис. 10, а и б.

Рис. 10. Диаграммы зависимости динамических перемещений сборного ротора от угловой скорости вращения: передняя опора (а) и задняя опора (б)

При анализе приведенных графиков сложно выявляется прохождение низшей критической скорости, что связано с низкой частотой вращения ротора при этом, и, соответственно, с низкими динамическими нагрузками. Явление резонанса выявляется при снижении углового ускорения раскрутки и «затягивании» прохождения через эту критическую скорость. Полученные низшие критические скорости сборного ротора составляют приблизительно 15 об/с, 80 об/с и 250 об/с. Соответствующие собственные частоты колебаний монолитного аналога ротора, вычисленные с применением математического аппарата извлечения собственных значений матрицы [9, 15], составляют приблизительно 25 Гц, 60 Гц и 280 Гц (см. таблицу).

Критические скорости вращения моделей роторов

Тип конструкции ротора Критические скорости вращения ротора, об/с (Гц)

1-я 2-я 3-я

Сборный 15 80 250

Монолитный 25 60 280

Несмотря на небольшое число поверхностей сопряжения, прослеживаются отличия сборного ротора от его монолитного аналога, указанные в работе [7]. Третья критическая скорость вращения сборного ротора меньше третьей собственной частоты колебаний монолитного ротора и располагается ближе к рабочему диапазону. Увеличение количества поверхностей сопряжения приводит к дальнейшему снижению жесткости сборного ротора и большим различиям в сравнении с результатами расчета монолитного аналога.

ЗАКЛЮЧЕНИЕ

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

СПИСОК ЛИТЕРАТУРЫ

1. Колотников М.Е. Предельное состояние деталей и прогнозирование ресурса газотурбинных двигателей в условиях многокомпонентного нагружения / под ред. В.М. Чепкина. -Рыбинск: Изд-во РГАТА, 2003. - 136 с.

2. Скубачевский Г.С. Авиационные ГТД, конструкция и расчет деталей. - М.: Машиностроение, 1981. - 552 с.

3. Хронин Д.В. Конструкция и проектирование авиационных газотурбинных двигателей. - М.: Машиностроение, 1989. - 565 с.

4. Леонтьев М.К. Современные методы расчета динамических характеристик роторных систем. Nastran или Dynamics? // Двигатель. - 2004. - № 3 (33). - С. 14-16.

5. Хронин Д.В. Теория и расчет колебаний в двигателях летательных аппаратов. - М.: Машиностроение, 1970. - 412 с.

6. Тимошенко С.П., Янг С.Х., Уивер У. Колебания в инженерном деле. - М.: Машиностроение, 1985. - 472 с.

7. Пыхалов А.А. Контактная задача статического и динамического анализа сборных роторов турбомашин: дис. .. .д-ра техн. наук: 05.07.05. - М., 2006. - 428 с.

8. Пыхалов А.А., Милов А.Е. Контактная задача статического и динамического анализа сборных роторов турбомашин: монография. - Иркутск: Изд-во ИрГТУ, 2007. - 192 с.

9. Зенкевич О. Метод конечных элементов в технике. - М.: Мир, 1975. - 542 с.

10. Расчеты машиностроительных конструкций методом конечных элементов: справочник / В.И. Мяченков, В.П. Мальцев, В.П. Майборода, В.Б. Петров, А.Н. Фролов, С.П. Заякин, Г.И, Ольшанская, В.Б. Горлов, В.С. Бондарь, С.П. Горшков, С.С. Корольков, Ю.В. Жуков, А.В. Цвелих; под общ. ред. М.И. Мяченкова. - М.: Машиностроение, 1989. - 520 с.: ил.

11. Александров А.В., Потапов В.Д. Основы теории упругости и пластичности: учебник для строительных специальностей вузов. - М.: Высшая школа, 1990. - 400 с.: ил.

12. Тимошенко С.П., Гудьер Дж. Теория упругости. - М.: Машиностроение, 1975. - 500 с.

13. Бабаков И.М. Теория колебаний: учебное пособие. - 4-е изд., испр. - М.: Дрофа, 2004. - 591 с.

14. NX Nastran. Handbook of nonlinear analysis (solutions 106 and 129). - [S. l.]: Siemens Product Lifecycle Management Software, 2014. - 661 p.

15. Бате К., Вильсон Е. Численные методы анализа и метод конечных элементов: пер. с англ. - М.: Стройиздат, 1982. - 448 с.

16. Дудаев М.А. Матрица жесткости балки Тимошенко в конечноэлементном анализе динамического поведения роторных турбомашин // Вестник Иркутского государственного технического университета. - 2014. - № 6 (89). - С. 59-65.

Дудаев Михаил Алексеевич, аспирант кафедры «Сопротивление материалов и строительная механика» Иркутского национального исследовательского технического университета. Основное направление научных исследований - метод конечных элементов. Имеет 8 публикаций. E-mail: [email protected].

Пыхалов Анатолий Александрович, доктор технических наук, профессор, кафедра «Сопротивление материалов и строительная механика» Иркутского национального исследовательского технического университета. Основное направление научных исследований - метод конечных элементов в механике твердого деформируемого тела. Имеет более 120 публикаций. E-mail: [email protected].

The contact problem in the analysis of the dynamic behavior of modular turbomachine rotors *

M.A. DUDAEV1, A.A. PYKHALOV2

1 Irkutsk National Research Technical University, 83, Lermontova Street, Irkutsk, 664074, Russian Federation, postgraduate. E-mail: [email protected]

2 Irkutsk National Research Technical University, 83, Lermontova Street, Irkutsk, 664074, Russian Federation, D.Sc. (Eng.), professor. E-mail: [email protected]

The article is about the dynamics of modern modular turbomachine rotors in transient operation modes. The relevance of solving contact mechanics in the analysis of the dynamic behavior of modular rotors of turbomachines is shown. The stiffness parameter of the modular rotor as a dynamical system is shown as a complex value that is a rotation speed function. The finite elements method is used to build a mathematical model of the modular rotor dynamics. The main stiffness components and forces of the modular rotor, their relation with physical processes characteristic of modular rotor systems are described in the paper. A mathematical model of solving the contact mechanics problem is proposed and a conventional algorithm to solve it using the finite element method is also described. Rotor excitation by force and moment imbalances of disks that make up its structure as well as imbalance modeling using the finite element model are considered. To solve the fundamental motion equation of the system by the finite element method a modified Newmark numerical integration method with a reduced mass matrix, a force vector and an adaptive integration time step is proposed. In so doing the problem is solved relative to an unknown acceleration vector at the time integration step. A test modular rotor finite element model with elastic supports that have contact interfaces in the zone of the shaft flange and disk connection is described. This model can be used for debugging software algorithms developed on the basis of the mathematical model proposed. The test model calculation results given as three minor oscillation forms and charts of dynamic displacements of the rotor in elastic support zones are shown. A brief analysis of the calculation results as well as a comparison of critical speeds of the modular rotor vibration frequency and natural vibration frequency of its monolithic analog is also presented.

Keywords: modular rotor, contact problem, finite element method, mathematical model, critical speed of rotation, dynamic behavior, Newmark method, transient process, vibrations

DOI: 10.17212/1814-1196-2015-3-113-129

Recivied 30 April 2015.

REFERENCES

1. Kolotnikov M.E. Predel'noe sostoyanie detalei i prognozirovanie resursa gazoturbinnykh dvigatelei v usloviyakh mnogokomponentnogo nagruzheniya [Limit state of parts and life prediction of gas turbine engines in a multi-component loading]. Rybinsk, RGATA Publ., 2003. 136 p.

2. Skubachevskii G.S. Aviatsionnye GTD, konstruktsiya i raschet detalei [GTE, construction and calculation of details]. Moscow, Mashinostroenie Publ., 1981. 552 p.

3. Khronin D.V. Konstruktsiya i proektirovanie aviatsionnykh gazoturbinnykh dvigatelei [The construction and design of aircraft gas turbine engines]. Moscow, Mashinostroenie Publ., 1989. 565 p.

4. Leont'ev M.K. Sovremennye metody rascheta dinamicheskikh kharakteristik rotornykh sis-tem. Nastran ili Dynamics? [Modern methods of calculation of dynamic characteristics of rotor systems. Nastran or Dynamics?]. Dvigatel' - Engine, 2004, no. 3 (33), pp. 14-16.

5. Khronin D.V. Teoriya i raschet kolebanii v dvigatelyakh letatel'nykh apparatov [Theory and calculation of vibrations in engines of aircraft]. Moscow, Mashinostroenie Publ., 1970. 412 p.

6. Timoshenko S.P., Yang S.Kh., Uiver U. Kolebaniya v inzhenernom dele [Vibrations in engineering]. Moscow, Mashinostroenie Publ., 1985. 472 p.

7. Pyhalov A.A. Kontaktnaya zadacha staticheskogo i dinamicheskogo analiza sbornykh roto-rov turbomashin. Diss. dokt. tehn. nauk [The contact problem of static and dynamic analysis of modular rotors of turbomachines. Dr. eng. sci. diss.]. Moscow, 2006. 428 p.

8. Pyhalov A.A., Milov A.E. Kontaktnaya zadacha staticheskogo i dinamicheskogo analiza sbornykh rotorov turbomashin. Monografiya [The contact problem of static and dynamic analysis of modular rotors of turbomachines. Monograph]. Irkutsk, ISTU Publ., 2007. 192 p.

9. Zenkevich O. Metod konechnykh elementov v tekhnike [The finite element method in the technique]. Moscow, Mir Publ., 1975. 542 p.

10. Myachenkov V.I., Mal'tsev V.P., Maiboroda V.P., Petrov V.B., Frolov A.N., Zayakin S.P., Ol'shanskaya G.I., Gorlov V.B., Bondar' V.S., Gorshkov S.P., Korol'kov S.S., Zhukov Yu.V., Tsvelikh A.V. Raschety mashinostroitel'nykh konstruktsii metodom konechnykh elementov. Spravoch-nik [Calculations of engineering structures using finite element method. A Handbook]. Moscow, Mashinostroenie Publ., 1989. 520 p.

11. Aleksandrov A.V., Potapov V.D. Osnovy teorii uprugosti i plastichnosti [Fundamentals of the theory of elasticity and plasticity]. Moscow, Vysshaya shkola Publ., 1990. 400 p.

12. Timoshenko S.P., Gud'er Dzh. Teoriya uprugosti [The theory of elasticity]. Moscow, Mashinostroenie Publ., 1975. 500 p.

13. Babakov I.M. Teoriya kolebanii [Theory of vibrations]. 4th ed., rev. Moscow, Drofa Publ., 2004. 591 p.

14. NX Nastran. Handbook of nonlinear analysis (solutions 106 and 129). Siemens Product Lifecycle Management Software, 2014. 661 p.

15. Bathe K.-J., Wilson E.L. Numerical methods in finite element analysis. Englewood Cliffs, New Jersey, Prentice-Hall, 1976. 528 p. (Russ. ed.: Bate K., Vil'son E. Chislennye metody analiza i metod konechnykh elementov. Translated from English. Moscow, Stroiizdat Publ., 1982. 448 p.).

16. Dudaev M.A. Matritsa zhestkosti balki Timoshenko v konechnoelementnom analize dina-micheskogo povedeniya rotornykh turbomashin [Timoshenko beam stiffness matrix in finite element analysis of turbomachine dynamic behavior]. Vestnik Irkutskogo gosudarstvennogo tekhnicheskogo universiteta - Bulletin of Irkutsk State Technical University, 2014, no. 6 (89), pp. 59-65.

ISSN 1814-1196, http://journals.nstu.ru/vestnik Science Bulletin of the NSTU Vol. 60, No.3, 2015, pp. 113-129