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

No 4

2006

534Л: 53 2.516.5+621.822.5

ДИНАМИКА РАЗГОНА ЖЕСТКОГО РОТОРА НА ПОДШИПНИКАХ

ЖИДКОСТНОГО ТРЕНИЯ

Канд. техн. наук А.О. ПУГАЧЕВ, д-р техн. наук проф. Л.А. САВИН,

канд. техн. наук, доц. О.В. СОЛОМИН

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

This paper describes the transient vibration of a rotor on flu id-film bearings. The model consists of a simple unbalanced rigid disc-shaft rotor and a hydrostatodynamic bearing. The rotor speed is time dependent to simulate startup processes. The system of motion equations includes the momentum equation. The driving moment is described by a simple expression. The fluid-film force is defined by integration of 2D Reynolds and the energy equations. The results of numerical calculations are presented in graphical form and discussed.

В конструкциях роторных систем высокоскоростных турбомащин различного назначения находят широкое применение опоры жидкостного трения [1—5]. Одной из важных задач, возникающих при проектировании и расчете системы «ротор — подшипники жидкостного трения», является обеспечение работоспособности агрегата на переходных режимах и, в частности, на режиме разгона, В первую очередь это относится к изучению проблемы нестационарного прохождения критических частот [6—9].

Рассмотрение этой задачи представляет значительный интерес для решения проблемы активного управления движением ротора [8] и реализации системы активной балансировки [6, 7]. Следует заметить, что отдельной задачей является выбор мощности привода роторной системы и/или закона управления вращающим моментом при использовании привода ограниченной мощности для разгона несбалансированного ротора до заданной скорости в необходимое время [7, 8] (во избежание проявления эффекта Зом-мерфельда). Кроме этого, определение динамических характеристик ротора при проходе критических частот позволяет оценить наибольшие напряжения в материале ротора, что важно для выполнения прочностного расчета роторной системы [9].

Характер развития динамических процессов в высокоскоростном малонагружен-ном роторе существенно зависит от закона разгона [10] и может стать причиной недопустимого роста амплитуд колебаний. Отметим, что процесс разгона вообще является для опор жидкостного трения, в особенности, гидродинамических, наиболее критическим режимом, поскольку сопровождается обкаткой или столкновениями цапфы вала и втулки подшипника и ведет к интенсивному износу последней. В ряде случаев, например, в кислородных турбонасосах ракетных двигателей [11], повышенная вибрация на переходных режимах и контакт цапфы и втулки могут привести к катастрофическим последствиям. Таким образом, представляется актуальным разработка модели и методов численного расчета динамики разгона системы «ротор— подшипники жидкостного трения».

№4

2006

Примем к рассмотрению ротор на двух одинаковых опорах, состоящий из симметричного невесомого жесткого вала и насаженного на него диска массой 2т и моментом инерции 27 (рис. 1). Отметим, что такая модель достаточно широко применяется в задачах динамики роторных систем [1—3, 7, 8]. Эксцентриситет положения центра масс диска (дисбаланс) равен Д. Не умаляя общности рассуждений, в качестве опоры рассмотрим полноохватный гладкий гидродинамический подшипник (ГДП) (рис. 2, а). Включе-

и

ние в расчетную схему других типов опор отразится только на расчете реакции смазочного слоя подшипников [12]. Ротор движется в плоскости радиального зазора опор с переменной угловой скоростью ш под действием: гидродинамических реакций опор /?, сил трения смазочного слоя моментов сил трения Мр крутящего момента 2А/К(сю) и веса 2В случае контакта цапфы и подшипника в уравнения движения роторной системы необходимо включить силовые факторы контактного взаимодействия: нормальную реакцию И* силу трения Т и момент трения ТВ!2 (рис. 2, б).

2т; 23

центр подшипника

центр масс ?/иДог

центр цапфы

Рис, 1. Динамическая модель жесткого ротора на опорах жидкостного трения

Процесс разгона ротора на ГДП включает три основных этапа:

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

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

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

Рассмотрим разгон ротора в ГДП после отрыва цапфы от втулки до достижения номинальной угловой скорости (третий этап). Используем уравнения Лагранжа [13], которые для обобщенных координат (X, К ф) примут вид

(1 дК дК

№4

2006

а)

б)

Рис, 2. Схема расчета ГДП: а — жидкостное трение; б — граничное трение

Кинетическая энергия К ротора может быть определена в виде [13]

2 т К; [и + 2тЛ2 )со2

к=-^ + -1—

2

2

(2)

Квадрат скорости центра масс определяется выражением (рис. 1)

V/ = X

2

+ У2 +2Дш(Усозф-Хз1пф) + Д2ш2.

(3)

Обобщенные силы (рис. 1) определяются следующими выражениями:

<2

X

2У/Х; {2у=2]Уу-2щ\ (2ф = 2М

к

2 М

/

2я^Дсо8ф,

(4)

где = /?х + ^ и1Уу = — составляющие полных реакций смазочного слоя; мно-

житель «2» в соотношениях (4) при и характеризует наличие двух одинаковых

опор. Подстановка выражений (2)—(4) в (1) дает

2 тХ = 2тД (ф2 соз ф + ф зт ф) + О,

2 т¥ = 2т А (ф2 з1п ф - ф соэ ф) + (2

(5)

(27 + 2шД2)ф = 2тД(х втф-Кссвф)*^

Используя первые два уравнения системы (5), исключаем линейные ускорения из третьего уравнения, что позволяет записать уравнения движения в виде

тХ = тД (ф2 соз ф + ф ем ф) + ^

тУ = тД (ф^тф-фсоБф)-!- -Уф = Ме + 1УхДзтф- М^Дсовф-М

Введем следующие безразмерные переменные и комплексы [12, 14]:

X

X

к

У

V

У

о

к

г

г

о

I

о

р0ОЬ

5

м

м,

Е .

Я

м

м

м

I

/

о

Ро° Ь

А

Л

тН

о

т

ЧРоОЬ

Л

]

3

'с Ро Р ~Ь

а

тЛ

Ц Р а®1'

О

те

р0ОЬ

М

М

(7)

р0В~Ь

номинальной

где г — время; = 2я/со.

угловой скорости ротора сон; Мн — номинальный крутящий момент привода ограничен-

ной мощности; к

о

радиальныи

<а; и и ь — соответ-давление подачи сма-

ственно диаметр и длина опорной поверхности подшипника; р0 —

зочного материала в подшипник.

Выражения (7) позволяют моделировать поведение роторной системы на основе критериев подобия и записать уравнения движения в безразмерном виде (штрих « ' » означает дифференцирование по безразмерному времени)

Ат X* = О. (ф/2 соэ ф + <р* вт ф) + А ¥" =

т

Л,ф' =

£2(ф'2 зтф- ф*со$

ММЁ + з1п <р-№гу соб ф)

(8)

М

/

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

контакте цапфы и втулки (рис. 2,6),

й

с

X

г^втос + Гсова); <2у =2(//со8а-Г8та); <2

с ф

ГА

(9)

цапфы

>

необходимо рассматривать неголономную систему, что обусловлено возможностью проскальзывания цапфы в точке контакта. В этом случае, строго говоря, для вывода уравнений движения необходимо воспользоваться уравнениями Аппеля [13]. Реакцию N и силу трения Т можно выразить через коэффициент трения для контактирующих поверхностей и коэффициенты жесткости и демпфирования при решении задачи теории упругости для внутреннего контакта двух цилиндров [15].

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

турбонасосов при условии достаточной мощности на турбине [16],

М

М0 (2 - со/0}

(10)

Чтобы определить реакции смазочного слоя Я, силы трения Г и момент Мг нужно знать, как распределяются давления р(х, г) по опорной поверхности подшипника' Для ГДП (рис. 2, а), смазываемого вязкой сжимаемой средой, поле давлений в условиях турбулентного течения находится из уравнения Рейнольдса [2, 14]

Э

/

э*

1г3 др

\

V

\лКх дх

+

в

У

дг

Г

Л3 др

\

ЦЛ\ дг

= 6

Э

Э.

(ргУ/г)-12рУ + 12/г-~

Эр

Э/'

(11)

2006

№4

где ц и р — вязкость и плотность смазочного материала; х иг — координаты по опорной

поверхности соответственно в окружном и осевом направлениях; Кх и енты турбулентности потока [2, 14]; Н(х, г) — функция радиального заз-скорости точек на поверхности цапфы. Выражения для определения зазор ростей и и V имеют вид (а = 2х/О)

эффици

- Хвта + Усова; и = а)Д+ Хсоза + УБта; У = Хеш а-У сое а.

(12)

В ряде практически важных случаев требуется учет тепловых явлений в смазочном слое [2, 7, 12]. В [14] предлагается использовать уравнение баланса энергий, записанное относительно энтальпии 1 (ср — теплоемкость при постоянном давлении смазочного материала; Т—температура):

Р

Э/ др дТ

--—+ с* —

др дг р дt

г

+ р

и

/Г др

\

\

2 12 \хКг дх

л V

д! др дТ --—+ с -—

др дх р дх

р/г2 др

12 цДГ Эг

д! др дТ

--£1 + с-

др дг дг

" I ~ I г- ЦЛ „.

дг 2 дх /г ^ А

(13)

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

Решение (11) и (13) дает распределение давлений по опорной поверхности подшипника. Интегрирование поля давлений р(х9 г) позволяет определить силовые факторы подшипника (т - напряжение сдвига в смазочной пленке)

пй I

710 Ь

Л

х

р$тас1гсЬс\ Яу = | ^ рсо$ас1г<1х\

/г Эр |хи X —---н

о о

по

о о

пО

К

х

11 тсо8 сикс!х\ ^ = -11теш ас1гс1х\ М

/

о о

о о

2 дх

В 2

/г

пй I

(14)

О о

Система дифференциальных уравнений (8), описывающая движение ротора, для

и

численного интегрирования приводится к системе уравнении первого порядка вида: у' = у Для ее решения предлагается использовать метод прогноза-коррекции Адам-

са—Башфорта—Моултона с адаптивным выбором шага интегрирования (и — начальные приближения на (п + 1)-ом шаге) [17]:

прогноз

у(0!

УП + 1

у + —(55/, л, 24 V

59/,1.1+37/я.2-9/„_з),

А/

коррекция уи+1 = Л + Т7 (9/),(°) +19/„ - 5/„_, - /„_,).

№4

2006

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

/

Входные величины: /г, с/, V, рк'\ Цикл 1: пока не выполнены условия сходимости

Расчет: р, ц, т, с

р

Решение уравнения баланса энергии: Г Цикл 2: пока не выполнены условия сходимости

р

Расчет: р, ц, т, с Расчет: кх, кг Решение уравнения Рейнольдса: р

Цикл 2; конец Цикл 1: конец

Рис. 3. Алгоритм итерационного определения поля давлений в ГДП

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

В соответствии с МКЭ представим функции давлений и энтальпий в виде

м

м

РИХР„Д„; I = /,„ Nln,

/И= 1

(16)

где /г , I

^ г ш5 m

давления и энтальпии в узлах сетки; М— число узлов сетки; N

U

ные базисные функции, которые для г-го узла имеют вид

N, =а(. +

(17)

ß

что N. равна единице в /-ом узле и обращается в нуль в других узлах.

Запишем (11) в безразмерном виде (здесь и далее в расчете давлений и энтальпий

используются уже безразмерные параметры р, /, /г, л, г [12, 14])

лд2р др д2 р др А—+ S —-t-C —^т + Е—

дх

дх

dz

dz

F

(18)

л r>

А = —В

h

Зр Эй ,3

+ h

' Р Л

дх

дх

\

11КХ

J-

D1 /г3р L2 ¡lKz

С

Е

D

L

h

Зр Эй , Э

+ й

Р Л

F

6D)x(

Pohl

р и

дh dx

\хК_ dz

+4(pi/)+

Эг

V

J-

Щ

>

2D

К

р V + h

\

Эр dt

\

JA

2006

Аналогично запишем уравнение энергий (13) в безразмерном виде [12, 14]

.31 . Э/ , Э/ _

Д — + А„-—А, —-С,

дt

Л Эх

Л

1

к

г

Р; = -^-р

Ро

V

1 и У^рц Н2 др

Г0 2 \10П \2\лКх дх

\

\ А.

/

^оУо Р^2 др

12^, Эг '

В

!

Ро др . в _ РоК У др

РоА/о &

Яро'сЛ 2 Эх

; 5

гу2^,.

РсЛ'

2 О

1г

; В — В1 + Вх + 5,;

С

л дТ л дТ л дТ 0

Э7

л- р

Эх

г р

Эг

(19)

Используя слабую формулировку метода взвешенных невязок по Галеркину [18] и аппроксимацию давлений в форме (16), сводим решение уравнения Рейнольдса (18) для треугольного элемента с вершинами (/, у, к) к решению системы линейных алгебраических уравнений вида (здесь р — вектор давлений [р., р., рк])

Кр = /.

(20)

Матрица К и вектор правой части для конкретного конечного элемента определя-

ются следующими соотношениями (А — площадь КЭ):

к = А

пт

/

А В В +5 Ь—С у V +Е У

пт * тп * т* п J

N

т

3

тп

з

\

У

А

(21)

Аналогично решение уравнение энергий (19) сводится к решению системы линей

ных алгебраических уравнений вида

г т

М — + К1 = /

йг

В данном случае элементы матриц Ми К я вектора/можно записать в виде

(22)

т

тп

А

I д2

тп

3

к

тп

д

3

("А-нл^л ^тпУт ) >

с

тт

А

3

(23)

Выражения (21) и (23) представляют собой соответствующие матрицы для конкретного КЭ. Сборка глобальных матриц для системы треугольных элементов, как правило, не представляет значительных сложностей. Процесс сборки подробно описан в [18]. Расчет по соотношениям (20)—(23) позволяет рассчитать поле давлений, которое после интегрирования дает реакции смазочного слоя (14).

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

АУМАНА

рр, ГИСТЕКА

№4

2006

что весьма информативным графическим представлением динамики переходных процессов в роторных системах, особенно для неустойчивых динамических процессов, является построение трехмерных изображений траектории движения центра цапфы в координатах X - У — со (рис. 4, б и 5, б) или X - У-г [20]. В этом случае достаточно наглядными становятся такие явления как прямая и обратная прецессия, контактные взаимодействия.

> о -

-0,2

-0.4 -

-0.6 -

-0.8 -

а)

б)

Тип Материал И, ь /го ро Р\ То т. Д и),1 мп

гдп Но 0.05 25 ■ Ю-0 0,75 ■ 10й т 0,25 • 10й 20 20 10 1000 5

Рис. 4. Разгон ротора на ГДП — орбитальное движение

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

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

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

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

№4

2006

>- О

X

О)

а)

б)

Тип Материал Б, Ь //о Р0 Р\ То т Д ........ М.,

ГДП Н) 4* 0,05 25 • 10-в 0,75 • 106 0,25 • 10й 20 15 10 0 2000 0,5

Рис. 5. Разгон ротора на ГДП — обкатка

3 0.5

ио 0.05 0,1 0.15 0.2 0.25 0.3 0.35 * 0.4 0.45 0.5

Тип Материал В 1 Ло Ро ■ ■ .'АР1 ■'.,' т '..'А *>

ГДП НгО '40 38 50 • 10~° »1 0,3- 10й 0,1 • 10" 297: 5: 5 - КГ6 700

Рис. 6. Эффект Зоммерфельда при разгоне ротора на ГДП

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

№4

2006

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

1. Давыдов А.Б.,Кобулашвили А. ILL, Ш е р ст ю к А. Н. Расчет и конструирование турбодетандеров.

— М.: Машиностроение, 1987. — 230 с,

2. А р т е м е н к о Н. П. и др. Гидростатические опоры роторов быстроходных машин, — Харьков: Основа,

1992.— 198 с.

3. Равикович Ю. А. Конструкции и проектирование подшипников скольжения агрегатов двигателей

летательных аппаратов. — М.: Изд-во МАИ, 1995. — 58 с.

4. Максимов В.А.,Баткис ПС, Трибология подшипников и уплотнений жидкостного трения высокоскоростных турбомашин. — Казань: ФЭН, 1998. — 430 с.

5. Геращенко Б. И. Динамика закритических роторов лопаточных машин. М.: Компания Спутник*, 2000.

— 250 с.

6. Zhou S., Shi J. The analytical imbalance response of Jeffcott rotor during acceleration // Journal of manufacturing science and engineering, 2001. — May. — Vol. 23. — P. 299—302.

7. Yamamoto Т., I s h i d a Y. Linear and nonlinear rotordynamics. A modern treatment with applications. — New York, John Willey&Sons, 2001. — 326 p.

8. Кельзон А. С,, M а л и н и н Л. М. Управление колебаниями роторов. — СПб.: Политехника, 1992. —

120 с.

9. ДиментбергФ. М. Изгибные колебания вращающихся валов. — М.: АН СССР, 1959. — 248 с. Ю.Шульженко Н. Г., Воробьев Ю.С. Численный анализ колебаний системы турбоагрегат —

фундамент. - Киев: Наукова думка, 1991. — 232 с. П.Шерстянников В.А. Исследование динамики роторов ТЫ А ЖРД // Двигатель. — 2002. — № 6. — С. 18—21.

12. С а в и н JI. А., С о л о м и н О. В. Динамика жесткого ротора на подшипниках скольжения, смазываемых криогенной жидкостью // Известия вузов. Машиностроение, 2004. —№4. — С. 27—38.

13.Уиттекер Э.Т. Аналитическая динамика, - Ижевск: Издательский дом «Удмуртский университет», 1999.

— 588 с.

14. Савин Л. А.,Соломин О.В. Расчет подшипников скольжения, работающих в условиях двухфазного состояния смазочного материала // Известия вузов. Машиностроение, 2004. — №2. — С. 36—42.

15. Д ж о н с о н К. Механика контактного взаимодействия. — М.: Мир, 1989. —-510 с.

16. Рудис М. А., Сафонов A.B. Анализ динамических характеристик роторов турбонасосных агрегатов // Авиакосмические технологии: Труды III международной научно-технической конференции. — Воронеж: Изд-во ВГТУ, 2002. — С. 147—152.

17. М эт ыоз Д. Г., Ф и н к К. Д. Численные методы. Использование MATLAB. — М.: Издательский дом «Вильяме», 2001. — 720 с.

18. Зенкевич О,, М о р г а н К. Конечные элементы и аппроксимация. — М.: Мир, 1986. — 318 с.

19. Программа расчета характеристик подшипников скольжения с криогенной смазкой («Подшипник-Крио-ген») /Савин Л.А., Соломин О.В., Пугачев А.О, и др. Свидетельство об официальной регистрации программы для ЭВМ № 2000610593.

20. С а г а Ъ е 11 i S., D е 1 р г е t е С., G е n t a G. Orbital tubes // Meccanica. — 1997. — № 32. — R 485—492.