Study of vibration sensitivity of a rotor of stationary gas-turbine unit on slide bearings to change of main settings of the system

Full Text

Введение Параметрическое исследование динамики ротора газовой турбины проводится на различных этапах её жизненного цикла - конструирование, сборка и наладка, эксплуатация и обслуживание. На этапе разработки перед инженером стоит задача в проверке пригодности выбранных технических решений. С другой стороны, во время эксплуатации интерес представляет анализ данных, полученных от системы мониторинга динамики ротора, и их интерпретация. В настоящей работе, в рамках создания модели динамики ротора, с применением стандартного программного комплекса МКЭ была разработана балочная модель ротора газовой турбины, с помощью которой были определены собственные частоты ротора на жестких опорах. На втором шаге исследования была разработана нелинейная модель ротора с учетом жесткостей нелинейных опор с подшипниками скольжения и протестирована с использованием модели ротора в стандартном программном комплексе МКЭ. В дальнейшем нелинейная модель ротора была применена для прямого интегрирования уравнений движения системы "ротор-подшипники" с использованием схемы интегрирования Ньюмарка 1. Модель опоры ротора с подшипниками на жидкой смазке Разработка модели опоры стационарной газотурбинной установки требует учета вклада в суммарную жесткость и демпфирование опоры ряда элементов, а именно жесткости корпуса установки одновременно с характеристиками жесткости и демпфирования смазки в подшипнике скольжения (рисунок 1). Характеристики жесткости и демпфирования смазки в подшипнике вычисляются с учетом деформаций рабочих поверхностей подшипника и вала. Для разработки модели опоры применена гипотеза о последовательном включении в единую расчетную схему упругих и демпфирующих элементов. Матрица жесткости специального конечного элемента опоры имеет вид: , (1) где: [Kb] ‑ матрица жесткости слоя смазки, учитывающая деформации поверхностей скольжения в подшипнике; [Kc] ‑ матрица жесткости корпуса опоры. Элементы матрицы демпфирования конечного элемента опоры [Csup] определяются с учетом компонент матрицы [Cb]. Рисунок 1. Подшипник скольжения: 1 – шейка вала; 2 – обойма подшипника Коэффициенты матрицы [Kc] определяются при расчете конструкции корпуса опоры. Расчет жесткости корпуса опоры может быть выполнен методом конечных элементов. Степень детализации при создании конечно-элементной модели зависит от особенностей конструкции. Для газотурбинных установок, корпуса которых часто представляют собой тонкостенные конструкции, при вычислении жесткостей опор необходимо моделировать весь корпус целиком. Для определения коэффициентов матриц [Kb] и [Cb] необходимо решать многодисциплинарную проблему упругогидродинамического контакта, включающую в себя задачу расчета деформаций рабочих поверхностей вала и подшипника и задачу течения смазки в зазоре. С учетом допущений для параметров течения несжимаемой смазки в подшипнике скольжения уравнения Навье-Стокса и уравнение неразрывности для слоя смазки [1] могут быть преобразованы к уравнению Рейнольдса [2]: (2) где оси координат выбраны следующим образом: s = φ·R – на одной из поверхностей скольжения в направлении относительного движения; z – на поверхности скольжения перпендикулярно к относительному движению; p(s,z,t) ‑ давление в слое смазки; h(s,z,t) – толщина слоя смазки; η ‑ вязкость смазки; R – радиус шейки вала; l – длина подшипника; ω – угловая скорость вращения; e(φ) радиус прецессии с угловой скоростью Ω; pa ‑ давление окружающей среды. В общем случае функция толщины слоя смазки h(s,z,t) определяется начальной формой и изменением взаимного расположения поверхностей скольжения (смещения шейки вала и повороты колодок в подшипнике) – hgeom, а также деформациями поверхностей скольжения под действием гидродинамических давлений [4] ‑ hdef . (3) При расчете h(s,z,t) учитываются взаимные смещения поверхностей скольжения в подшипнике, как абсолютно жестких, и деформации поверхностей скольжения, вызванные давлением смазки p(s,z,t) в зазоре. После процедуры конечно-элементной дискретизации [3-5] система уравнений (2) сводится к системе конечно-элементных уравнений для определения давления [2]. Определение деформаций поверхностей скольжения подшипника под действием давления смазки требует, в общем случае, решения трехмерной задачи теории упругости, которое осуществляется при помощи метода конечных элементов. Для создания трехмерных моделей шейки вала и подшипника используются объемные 20-ти узловые конечные элементы [5]. Типовая конечно-элементная схема для шейки вала и подшипника с гладкими поверхностями представлена на рисунке 2. Рисунок 2. МКЭ модель и деформации поверхностей скольжения подшипника Решение задачи ведется последовательными приближениями. На первом шаге подшипник и шейка вала считаются жесткими (hdef = 0). В этом случае из решения (2) следует распределение давления в смазочном слое. При приложении этого давления на соответствующие поверхности в трехмерные модели подшипника и вала, определяется первое приближение для деформаций рабочих поверхностей и, соответственно, hdef в первом приближении. Корректируя в (3) толщину слоя смазки, из (2) определяем новое распределение давления на втором шаге и затем второе приближение для деформаций поверхностей и т.д. Итерационный процесс останавливается, когда изменение давления и перемещений поверхностей на двух последующих приближениях находятся в пределах требуемого допуска. 2. Модель динамики ротора Для проведения исследований разработана многодисциплинарная модель динамики ротора, включающая в себя модель вала, вращающихся дисков и опоры, в которой для вала использована балочная конечно-элементная модель с распределёнными инерционными параметрами и учётом инерции поворота сечений при изгибных деформациях ротора. Масса и моменты инерции дисков учитываются конечными элементами сосредоточенной массы в заданных сечениях вала. Модель опоры представлена специальным конечным элементом, коэффициенты которого определяются в ходе вычислений для текущего положения ротора и скорости его вращения. Модель ротора в неподвижной системе координат, описывается уравнениями МКЭ следующим образом , (4) где: [M] – матрица масс вала и вращающихся деталей; [C] – матрица, учитывающая влияние гироскопических моментов и демпфирования в опорах и уплотнениях; [M1] – часть матрицы [M] связанная с линейными перемещениями узлов; , и - векторы перемещений узлов ротора, их скоростей и ускорений соответственно; [KR] – матрица жесткости ротора; [KS] – матрица жесткости, учитывающая влияние опор и уплотнений; {F0} – вектор внешних сил, действующих на ротор; {E0} – вектор начальных дисбалансов поперечных сечений вала, которые определяются {ε0}. Матрицы массы [M] и жесткости [KR] не изменяются в процессе движения с постоянной скоростью вращения. В противоположность им, коэффициенты матриц [KS] и [C] зависят не только от скорости вращения ротора, но также от положения искривленной оси вала в опорах. Коэффициенты матрицы [KS] определяются свойствами течения смазки в подшипниках и уплотнениях, аэродинамическими силами в колёсах и жесткостью корпуса ГТУ. Гироскопические моменты, параметры течения смазки в уплотняющих элементах и трение в смазочной плёнке в подшипниках влияют на коэффициенты матрицы демпфирования [C]. Таким образом, [KS] и [C] являются матрицами общего вида с коэффициентами, которые нелинейно зависят от текущего положения оси вала, скорости вращения и параметров масла. 3. Расчет параметров модели ротора с нелинейными опорами На рисунке 3 представлена типовая конструкция валопровода стационарной ГТУ, состоящая из двух валов соединённых муфтой. Каждый вал опирается на два радиальных гидродинамических подшипника. При разработке модели валопровода созданы соответствующие системе (4) балочная модель ротора с нелинейными опорами и балочная модель ротора в стандартном программном комплексе МКЭ, предназначенная для тестирования балочной модели с нелинейными опорами. Эти модели использовались для исследований динамики ротора при изменении параметров ГТУ. П4 П3 П1 П2 Рисунок 3. Модель валопровода стационарной газотурбинной установки 3.1. Тестирование модели ротора Верификация модели ротора с податливыми нелинейными опорами проводилась для жёстких опор и для опор ротора с жёсткостью 8 МН/мм. Результаты сравнения собственных частот для двух моделей представлены в таблице 1. Частоты для обеих моделей соответствуют друг другу. Собственные формы ротора представлены на рисунке 4. Таблица 1 Сравнение моделей № Собственные частоты, Гц Жёсткие опоры Опоры с жёсткостью 8 МН/мм Стандартный комплекс МКЭ Балочная модель Разница % Стандартный комплекс МКЭ Балочная модель Разница % ft1 22,10 22,12 0,09 21,68 21,70 0,09 ft2 34,35 34,51 0,47 32,23 32,35 0,37 ft3 56,31 56,77 0,82 54,86 55,29 0,78 ft4 93,32 94,34 1,09 81,81 82,72 1,11 a) ft1 – ротор 2 б) ft2 – ротор 1 Рисунок 4. Собственные формы ротора 3.2 Характеристики подшипников скольжения Расчеты параметров подшипника скольжения проводились для подшипника с гладкими поверхностями. Распределение давления в подшипнике с деформируемыми поверхностями, определенными с использованием объемной модели подшипника, представлены на рисунке 5. Рисунок 5. Распределение давления в подшипнике с податливыми поверхностями Линия “q” на рисунке 5 соответствует месту, где толщина слоя смазки минимальна для подшипника с жесткими поверхностями. Давление смазки в подшипнике с деформируемыми поверхностями заходит за линию “q”. В данном случае, когда оси вала и подшипника остаются параллельными при смещении шейки вала, распределение давления симметрично относительно линии “s”, проходящей через середину подшипника. Если имеют место угловые перемещения шейки вала в подшипнике, то распределение давления в подшипнике становится несимметричным относительно линии “s” и помимо подъемной силы, также возникает момент давления смазки, создаваемый относительно плоскости симметрии подшипника (линия “s”). Как следует из представленных результатов, учет податливости поверхностей изменяет картину распределения давлений в подшипнике, что сказывается на величине и направлении подъемной силы и, соответственно, на характеристиках жесткости и демпфирования слоя смазки. На рисунке 6 приведены зависимости безразмерной подъемной силы от значения относительного эксцентриситета и зависимость угла, характеризующего отклонение подъемной силы от направления эксцентриситета. Рисунок 6. Статические характеристики подшипника скольжения подъемная сила и угол между направлениями подъемной силы и эксцентриситета В рассматриваемом случае поле давлений симметрично относительно плоскости “s”, момент давления равен нулю и матрицы жесткости [Kb] и демпфирования [Cb] имеют размерность 2х2. 4. Исследование динамики ротора на нелинейных опорах Для решения системы уравнений движения (4) применена схема интегрирования Ньюмарка, причем на каждом временном шаге методом последовательных приближений выполняется уточнение коэффициентов матриц жёсткости и демпфирования. 4.1 Влияние конструкционного демпфирования Исследование влияния конструкционного демпфирования на динамику валопровода выполнено для значений коэффициента конструкционного демпфирования b: 0.001, 0.002 и 0.02 при условии что , где - матрица конструкционной жёсткости системы. Расчеты проведены для частот вращения валопровода от 0 до 80 Гц с шагом в 5 Гц. Полные спектры и орбиты валопровода для конструкционного демпфирования 0.001, 0.002 и 0.02 представлены на рисунках 7-8. a) b = 0.001 б) b = 0.002 в) b = 0.02 Рисунок 7. Орбиты ротора в подшипнике П1 (60Гц) и спектр с различным конструкционным демпфированием и вязкостью смазки 0.04 МПа·с a) b = 0.001 б) b = 0.002 в) b = 0.02 Рисунок 8. Орбиты ротора в подшипнике П3 (60Гц) и спектр с различным конструкционным демпфированием и вязкостью смазки 0.04 МПа·с Для случая b = 0.001 в интервале от 20 Гц до 55 Гц имеют место значительные колебания на частотах, близких к ft2, 2·ft2 и 3·ft2. В отличие от b = 0.001, для b = 0.002 значительные колебания ротора на частотах, близких к ft2, 2·ft2 и 3·ft2, появляются только на небольшом интервале между первой и второй критической скоростями с амплитудами меньшими в три раза. На всех остальных частотах ротор колеблется с малыми амплитудами. Для случая b = 0.02 колебания ротора на частотах близких 1X (прямая прецессия) имеют место для скоростей вблизи второй критической скорости ротора 2 (ft3). Характер колебаний вала в П3 на рабочей частоте 60 Гц практически не отличается от колебаний в П1. Однако для П3 при b = 0.001 значительные колебания на частотах близких к ft1 присутствуют только в интервале от 20 Гц до 25 Гц. 4.2 Влияние вязкости масла Исследование влияния вязкости смазки в подшипниках скольжения на динамику валопровода проведено для вязкости h равной 0.02 МПа·с, 0.04 МПа·с и 0.06 МПа·с. Расчеты проведены для частот вращения валопровода от 0 до 80 Гц с шагом в 5 Гц. Полные спектры и орбиты валопровода для вязкостей смазки 0.02 МПа·с, 0.04 МПа·с и 0.06 МПа·с представлены на рисунках 9-10. a) h = 0.02 б) h = 0.04 в) h = 0.06 Рисунок 9. Орбиты ротора в подшипнике П1 (60 Гц) и спектр при различной вязкости масла и конструкционным демпфированием 0.002 a) h = 0.02 б) h = 0.04 в) h = 0.06 Рисунок 10. Орбиты ротора в подшипнике П3 (60 Гц) и спектр при различной вязкости масла и конструкционным демпфированием 0.002 При вязкости смазки 0.02 МПа·с наблюдаются существенные вибрации ротора 1 на частотах близких к ft2, 2·ft2 и 3·ft2 для частот вращения ротора в интервале от 30 до 65 Гц. Ротор 2 вибрирует в этом интервале с амплитудами в 10 раз меньшими чем амплитуды ротора на его первой критической скорости (ft1). При вязкости смазки 0.04 МПа·с присутствуют существенные колебания на частотах, близких ft2 и 2·ft2 для частот вращения в интервале около первой критической скорости ротора 1 (ft2). Амплитуды колебаний ротора 2 в 10 раз меньше чем аналогичные для ротора 1 для интервала, начинающегося с частоты вращения ротора, равной ft2. Для вязкости 0.06 МПа·с колебания роторов соответствуют синхронной прецессии (±1X) с малыми амплитудами одного порядка. 4.3 Чувствительность к изменению дисбаланса Исследование влияния начального дисбаланса на динамику валопровода проведено для значений дисбаланса: 5, 10 и 15 кг∙м для всех частот вращения. Проведены расчеты для вращения валопровода на частотах в диапазоне от 0 до 80 Гц с шагом 5 Гц. Орбиты движения валопровода в подшипниках и их спектры представлены на рисунках 11-12. При дисбалансе равном 10 кг∙м валопровод колеблется на частотах кратных частоте вращения ±1X, ±2X и ±3X, амплитуды колебаний ротора 1 достигают 150 мкм, а для ротора 2 составляют менее 30 мкм. При дисбалансе равном 15 кг∙м валопровод колеблется на частотах кратных частоте вращения ±1X, ±2X и ±3X. Кроме того возникают субгармонические колебания на частоте 9 Гц и кратных ей частотах. На спектрограмме можно увидеть серию пиков соответствующих этим субгармоническим частотам. a) 5 кг∙м б) 10 кг∙м в) 15 кг∙м Рисунок 11. Орбиты ротора в подшипнике П1 (60 Гц) и спектр при различных дисбалансах a) 5 кг∙м б) 10 кг∙м в) 15 кг∙м Рисунок 12. Орбиты ротора в подшипнике П3 (60 Гц) и спектр при различных дисбалансах Как следует из сравнения амплитуд колебаний роторов в подшипниках П1 и П3 возбуждение первого ротора приводит лишь к малому по амплитуде возбуждению второго ротора. Выводы Разработана методика моделирования динамики ротора на нелинейных опорах. Модель ротора состоит из балочной модели вала с присоединёнными дисками и модели опорных подшипников скольжения с масляной смазкой, что позволяет вычислять собственные значения и формы, устойчивые и неустойчивые орбиты ротора для различных параметров модели. На базе модели динамики ротора исследована чувствительность валопровода стационарной ГТУ к изменению конструкционного демпфирования, вязкости смазки и величин дисбаланса роторов. Продемонстрировано влияние этих параметров на орбиты (амплитуды) и частоты колебаний ротора. Определены границы устойчивых и неустойчивых режимов работы в зависимости от параметров системы. Методика позволяет за счет изменения параметров системы, таких как конструкционное демпфирование, вязкость масла и дисбаланс, выявлять нежелательные режимы работы ротора на всех этапах жизненного цикла стационарной газотурбинной установки и определять методы ограничения амплитуд колебаний ротора для обеспечения устойчивости его вращения.

About the authors

Y. M Temis

Central Institute of Aviation Motors

Email: tejoum@ciam.ru
Dr.Eng., Prof.; +7 (495) 361-64-82

M. Y Temis

Central Institute of Aviation Motors

Email: tejoum@ciam.ru
Ph.D.; +7 (495) 361-64-82

A. M Egorov

Central Institute of Aviation Motors

Email: tejoum@ciam.ru
+7 (495) 361-64-82

References

Supplementary files