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

МАТЕМАТИЧЕСКАЯ МОДЕЛЬ ГИБКОГО РОТОРА НА ОСНОВЕ ОБОБЩЕННЫХ ЛАГРАНЖЕВЫХ КООРДИНАТ

В. Г. Быков1, А. Е. Мельников2

1. С.-Петербургский государственный университет, канд. физ.-мат. наук, доцент, vgbykov@mail.ru

2. С.-Петербургский государственный университет, аспирант, m02mae@mail.ru

Исследования динамики роторов насчитывают более чем 140-летнюю историю, о чем свидетельствует работа У. Рэнкина [1]. Первые математические модели вращающегося вала были созданы Фепплем [2] и Джеффкоттом [3]. Значительный вклад в теорию гибких роторных систем по таким актуальным проблемам, как расчет критических скоростей, балансировка, устойчивость движения внесли А. Стодола, Р. Граммель, А. Н. Крылов, Ф. М. Диментберг и др. Фундаментальными работами, отражающими современное состояние теории жестких и гибких роторов, являются монографии Генты [5] и Ямамото [6].

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

1. Описание модели и вывод уравнений. Рассмотрим ротор в виде тонкого диска, несимметрично насаженного на упругий вал круглого сечения с равномерно распределенной массой. Вал вращается в шарнирных опорах, одна из которых является подвижной (рис. 1). Диск считаем абсолютно твердым однородным телом массы т, жестко скрепленным с валом в точке 0\ —центре диска. Неуравновешенность ротора зададим с помощью точечного дисбаланса массы тр, закрепленного на диске в точке В Обозначим через е расстояние от центра диска до точки дисбаланса.

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

где Хк (і) и ук (і) —величины максимального прогиба вала, а I —длина вала.

Для описания движения диска введем подвижную систему координат Оіхіуїхі, оси которой параллельны осям инерциальной системы Охух, и систему координат О1Х2У2Х2,

© В. Г. Быков, А. Е. Мельников, 2010

(1)

к = 1

к=1

Рис. 1. Диск на гибком валу.

Рис. 2. Углы поворота диска.

жестко связанную с диском ротора. Ось O\Z2 направим вдоль нормали диска, а ось Ox —так, чтобы она проходила через точку D (рис. 2).

Обозначим через n вектор нормали к плоскости диска. Пусть р — угол между проекцией n на плоскость OiXizi и осью Oizi, а ф —угол образованный проекцией n на плоскость Oiyizi и осью Oizi (рис. 2,a). Переход от системы Oixiyizi к системе OiX2y2Z2 можно осуществить тремя поворотами: на угол р вокруг оси Oiyi, на угол ф вокруг OiX2 и на угол х вокруг оси Oiz2 (рис. 2,b). Угол х характеризует собственное вращение диска вокруг его нормали. Считая углы р и ф малыми, запишем матрицу перехода от системы Oixiyizi к системе OiX2y2z2:

{cos х sin х Р cos х + Ф sin X \

— sin х cos х —р sinх + Ф cos х / • (2)

-р —ф 1 )

В силу симметрии оси системы OiX2y2z2 совпадают с главными осями инерции диска.

Угловая скорость диска зависит от угловой скорости вращения вала и от деформации вала. Вектор-столбец проекций вектора угловой скорости диска в системе координат OiX2y2z2 имеет вид

0 = Р{-ф,р,0}Т + {0, 0,х}Т = jpsin х — /cos х, pcos х + /sin х,Фр - РФ + х} • (3)

Предполагая, что вал мало отклоняется от недеформированного положения, будем считать, что проекция центра диска на ось Ох будет оставаться неизменной. Обозначим через го1 = {X(хо,і),У(хо,і),хо}т вектор столбец проекций радиуса-вектора точки

Ох в неподвижной системе координат, а через г'в = {е, 0,0}т — радиус-вектор точки Б в системе координат О1Ж2у2X2. Тогда положение точки Б в неподвижной системе координат найдем по формуле

г Б = ГОї + РТ г'в

X (хо,Ь) + е 008 х У (хо,Ь) + е 8Іп х хо + е(р 008 х + Ф 8ІП х)

(4)

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

( дХ(х,і) \ дх

дУ (х,і)

дх

1

(5)

Так как углы р и ф считаются малыми, мы можем приближенно положить

р = arctg |

ф = аг^

дХ(х, і)

дх ' дУ(х,і)

дХ(х, і)

дх

дх дУ (х, і)

дх

™ кп кпхо

= Е -гхи^) сои,

2=20 к = 1 1 1

™ кп кпхо

= Е -гУкУ) сов——.

2=20 к = 1 1 1

(6)

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

Далее мы будем рассматривать конечное число п собственных форм колебаний. Введем следующие обозначения:

а = {ах, а2,..., ап}Т , Ь = {61, Ь2,..., Ьп}Т ,

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

кпхо кп кпхо

ак = вт —-—, 6,% = —сое—.

Для вывода уравнений движения ротора воспользуемся методом Лагранжа. Кинетическая энергия системы складывается из кинетических энергий вала, дисбаланса и диска:

1 1 1

Т = 2тя X 0^(*)2 + Ук{1?) + + 2 (т^О! + вТ(70)) . (7)

и = 1

где = рБ1/2 — приведенная масса вала; 3 = diag{3%, 3%, Зр} — матрица центрального

тензора инерции диска; Зр, 3% —соответственно, полярный и транверсальный моменты инерции диска.

п

г = 2п

2=20

2 = 20

2 = 20

Перейдем к матричной комплексной переменной V = х + iy, где х и у — вектор-столбцы, составленные из п соответствующих функций хь (Ь) и уь(Ь). Обозначив через V вектор комплексно сопряженный к V, кинетическую энергию системы с точностью до малых второго порядка представим в виде

1 • т 1

Т = -V (теЕ„ + + (то + тв)А) V + ~(1р + тпе2)х2-

— ,]рх\т\уТ Ву] + трех 1т[\гТе_*х]а, (8)

где Еп —единичная матрица п х п; А = а • аТ; В = Ь • ЬТ.

Так как на вал не действуют внешние потенциальные силы, потенциальная энергия ротора будет состоять только из энергии изгиба вала:

П = -шдутГ22у, (9)

где П2 = diag{иf,и'^ • • -,и2п}; <ш\ = k4п4EJ/(рБl4) — собственные частоты поперечных колебаний вала; р, Б, I, Е.1 — плотность, площадь поперечного сечения, длина и изгиб-ная жесткость вала соответственно.

Для учета вязкого трения введем в рассмотрение диссипативную функцию Релея

1 Т 1

я = 2С^ * + 2С^2’

где с и сх — коэффициенты сопротивления при изгибе и вращении вала.

С учетом выражений (8), (9) и (10) запишем уравнения Лагранжа 2-го рода:

{(швЕп + JtB + (ш + шп)А) V + (сЕ„ + 2iJpxВ^+

+ (швП2 + iJpxB)v = шве(х2 — ^)егха, (11)

(+ тде2)х + тоде 1т[уТе_гх]а + Мр 1т[у Ву] + сх^: = Мхт где Мх — внешний вращающий момент, приложенный к валу.

2. Стационарные режимы вращения ротора. Предположим, что ротор вращается с постоянной угловой скоростью х = и. Перейдем к новой векторной переменной w относительно системы координат Охзузх, вращающейся с угловой скоростью и вокруг оси Ох. Переменные V и w связаны соотношениями

V = , V = (У + шу)егш1^, V = (У + 2^У — и2у)егш*' • (12)

Подставляя соотношения (12) в уравнения (11), получаем

(шв Еп + JtB + (ш + шв)А) (У + 2^У — и2у) + (сЕп + 2iJpшB)(лv + ^у)+

+ шв П2у = шв еи2а. (13)

Введя новые переменные х, = х/1, У = у/е и т = и\Ь, представим уравнение (13) в безразмерном виде:

^ ., / С2У СУ 2 _ \ . _ ~ СУ _ N ~о 2

(Е„ +7(В + тА) I 2 -\-2tuj—-оглу I +(сЕ„ + 2* 1ршв) I ——|-*иш I + П лу = шци а,

Т Т Т (14)

J n'2Jt J n2Jp ... , ..o

Jt = -77--------, Jp = -77——, m = -, mjj

n2Jp m + mo mo

;9 , -р=ю---------------------------> m= -> mD = -.

l2ms l2ms ms ms

П2 ~ l2

ms wi wi w2 ’ n2

Далее символ «~» писать не будем.

Рассмотрим стационарные режимы вращения ротора, при которых w = w* = const. Полагая в (14) значения всех производных от w равными нулю, получаем линейную алгебраическую систему относительно n компонент вектора w*:

(П + w(ic — w)En — w (mA + (Jt + 2Jp)B)) w* = "mow a. (15)

Обозначим n x n-матрицу левой части системы (15) через R. Заметим, что эту матрицу можно записать в виде

R = Л + U • VT, (16)

где Л = П2 + w(ci — w)En —диагональная матрица размера n x n, а U = {a, b} и V= —w2 {ma, (Jt + 2Jp)b} —матрицы размера n x 2.

Будем искать матрицу, обратную к R, в виде

R-1 = Л-1 + Л-1 • U • D • VT • Л-1, (17)

где D —неизвестная матрица размера 2 x 2. Для того чтобы R-1 • R = R • R-1 = En, необходимо и достаточно выполнения условия

E2 + D +(VT • Л-1 • U) • D = 0. (18)

Условие (18) представляет собой матричное уравнение порядка 2x 2 относительно неизвестной матрицы D. Решая это уравнение и предполагая, что матрица E2 + VT • Л-1 • U не вырождена, находим

D = —(E2 + VT • Л-1 • U)-1.

Таким образом, вычисление обратной матрицы системы (15) сводится к вычислению диагональной матрицы Л-1 и определению обратной матрицы размера 2 x 2:

R-1 = Л-1 — Л-1 • U • (E2 + VT • Л-1 • U)-1 • VT • Л-1. (19)

Выражение (19) позволяет разрешить систему (15) относительно вектора w*:

w тр((» + 2Л^~5и)Л"‘а + 5ігЛ~1Ь) (20)

где

Sii = aT Л-Ч S2i = Si2 = aT Л-іЬ, S22 = bT Л-іЬ.

Для нахождения амплитудно-частотных характеристик (АЧХ) ротора определим величину прогиба оси вала в точке крепления диска:

(A(uj))2 = + (^2yk(t)a^ = (w*Ta)(w*Ta), (21)

Т _ 2 (____________1 — #22 (^ + 21р)ш2____________ 2 А /00\

а-'твШ ии(1-522(Л + 27р)с;2)+522(74+27р)с;2 та; ^ ' (22)

Аналогично, чтобы найти суммарную амплитуду углового отклонения диска, определим наклон нормали к поверхности диска:

/ ^ \ 2 / ^, \ 2

(Ф(ш)) = xk(t)bkj + ^7 ,yfc(t)b^ = (w*Tb)(w* b), (23)

где

Th = ____________________________S12mDuj2____________________________

W* (1 — Sjpii)2) (1 — S'22(Jt + 2Jp)lv2) — S22m(Jt + 2Jp)uji " 1 j

Предполагая теперь, что n ^ ж, представим выражение для Sii в виде ряда

то . 2 I

с T*-1 sln knz0 /ое-\

Sii = а Л а= > —----------;-----(25)

11 ^ к4 + аш — и2 v '

k=i

Чтобы найти сумму ряда (25), воспользуемся формулой

■ТО coskz 1 n (cos(a(n — z)) ch(a(n — z))

^ k4 — a4 2a4 4a3 V sln an sh an

k=1

которая верна для любых a£CH0<z<2f [10]. Подставляя в (26) z = 2nzo и а = \]ю2 — сш, получаем

ТО . 2 7 ч ТО 1 ч ТО 1 л - 27

Esin knzo 1 1 1 1 — 2sin knzo

kA + cSlo — lo1 2 kA — a4 2 kA — a4

k=i k=i k=i

n , n /cos(an(1 — 2zn)) ch(an(1 — 2zo)) \

= -Г-з ctgaTT + cthcwr ^— + -------— • 27

8a3 8a3 \ sin an sh an J

Остальные Sj найдем, проделав аналогичные преобразования. В результате будем иметь

n (sin anzo. sh anzo

йп = —5- -----sin cml — zo)----------sh cml — zo)

4a3 \ sin an sh an

n cos anzo ch anzo

612 = o' —:-----sin дат(1 - zo)--------sh cml - z0)

4a2 \ sin an sh an

n cos anzo ch anzo

022 = —-------:-----cosa7r(l - zo) H-------cha7r(l - z0) ) .

4a sin an sh an

В случае zo = 1/2, что соответствует закреплению диска посередине вала, формула (22) принимает более простой вид:

T 2 / 8a 2 \ / \

w* а = mDu а7Т--^ - ш . (28)

^(tg—-th-) J

к = 0.5 = 0.3

Прибл. Точн. Прибл. Точн.

1 0.25 0.4156 0.1276 0.4219

2 6.25 6.5545 0.6944 1.0000

3 20.25 20.5717 3.1888 3.3448

4 42.25 42.5781 10.3316 10.4570

5 72.25 72.5815 17.3611 17.6419

Уравнение для определения критических частот Vk найдем, приравняв в (22) знаменатель нулю и подставив а = в выражения для :

----------1~£1М + ЫрУк---------------_ ^2 (29)

Б11 (1 — Б22(Jt + 2‘1р)и‘2) + Б12(^ + 2Jp)v‘к

^Р)ик) ^ Б12\*ир)"к

Заметим, что если записать уравнение (29) в виде

( 2 — ^Ц4) ( / Т , 0 Т N 2 — $22^) — ^12;

\ш^к / \(Jt + 2'^р)^к /

то при достаточно больших Vк слагаемыми порядка 1/у\ можно пренебречь, а гиперболическая часть в выражениях для Б^ будет приближенно равна 0-5. Таким образом, для определения критических частот высокого порядка имеем простую приближенную формулу

со§((1 - 2г0)л/щ,-к) = вт(уед,

откуда находим

^ = (2 ±1(2_—4го))2 ’ ‘ = “.1.2.- (30)

3. Численное моделирование. При проведении расчетов были выбраны следующие значения параметров: I =1 м, 2шв = 3 кг, ш = 7 кг, шв =0.1 кг, К = 0.25 м, г = 0.01 м, е = 0.05 м, Е = 0.2 • 1010 кг/(м с2), с = 16.22 кг/сек.

В таблице представлены рассчитанные по формулам (29) и (30) точные и приближенные значения критических частот для случаев симметричного (хо = 0.5) и несимметричного (хо = 0.3) закрепления диска.

Из таблицы видно, что относительная погрешность приближенных значений, начиная с третьей частоты, не превышает 5%.

На рис. 3 показаны графики зависимости амплитуды отклонения центра диска от угловой скорости, рассчитанные по формуле (21) (сплошные линии), и аналогичные кривые для модели ротора с невесомым валом (пунктирные линии), приведенные в [11]. Рисунок 3^ соответствует симметричному закреплению диска, а рис. 3,Ь — несимметричному (хо = 0.31). На рис. 3,с для несимметричного ротора показана суммарная амплитуда углового отклонения диска. Для симметричного ротора данная амплитуда, определяемая формулами (23), (24), будет тождественно равна нулю, поскольку в силу симметрии у вала остаются только нечетные формы колебаний.

Анализ графиков показывает, что АЧХ для моделей ротора с весомым и невесомым валом практически совпадают в окрестности первой критической скорости, но в дальнейшем мы видим существенные количественные и качественные различия, связанные

А(ш)

А(ш)

Ф(ш)

с 0,10 0,08 0,06 0,04 0,02

0

1 2 3 4ш

Рис. 3. АЧХ с учетом и без учета массы вала.

1 2 3

Рис. 4. Переход через критические скорости.

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

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

На рис. 4 представлены амплитудные кривые при нестационарном прохождении через первые две критические скорости ротора с симметрично (рис. 4,a) и несимметрично (рис. 4,b) закрепленным диском. На рис. 4,c показана амплитуда углового отклонения диска в случае его несимметричного закрепления.

Графики на рис. 4 наглядно демонстрируют адекватность предложенной модели, а также обоснованность сделанных предположений о малости углов и отклонений.

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

Литература

1. Rankine W. J. McQ. On the Centrifugal Whirling of Shaft // The Engineer. Vol. 27. 1869. P. 249.

2. Foppl A. Das Problem det Laval’schen Turbinenwelle // Der Civilingenieur. Vol. 41. 1895. S. 333-342.

3. Jeffcott H. H. The Lateral Vibration of the Loaded Shafts in the Neighbourhood of a Whirling Speed // Phil. Mag. Vol. 6. N37. 1919. P. 303-314.

4. Диментберг Ф. М. Изгибные колебания вращающихся валов. М.: Изд. АН СССР. 1959. 248 с.

5. Genta G. Dynamics of Rotating Systems. Springer, 2005. 658 с.

6. Yamamoto T., Ishida Y. Linear and Nonlinear Rotordynamics: A Modern Treatment with Applications (Wiley Series in Nonlinear Science). Wiley-Interscience. 2001. 348 p.

7. Зегжда С. А., Юшков М. П. Применение уравнений Лагранжа первого рода при исследовании собственных колебаний вала с дисками // Механика твердого тела. 1999. №4. С. 31-35.

8. Зегжда С. А., Солтаханов Ш.Х., Юшков М. П. Уравнения движения неголономных систем и вариационные принципы механики. Новый класс задач управления / Под ред. проф. П. Е. Товстика. М.: ФИЗМАТЛИТ, 2005. 272 с.

9. Вернигор В. Н. Определение собственных частот и эквивалентных масс упругого тела по его динамической податливости // Вестн. Ленингр. ун-та. Сер. 1. 1990. Вып. 4 (№2). С. 35-42.

10. Градштейн И. С., Рыжик И. М. Таблицы интегралов, и сумм, рядов и произведений. М.: ФИЗМАТЛИТ, 1963. 1108 с.

11. Быков В. Г. Стационарные режимы движения неуравновешенного ротора с автобалан-сировочным механизмом // Вестн. C.-Петерб ун-та. Сер. 1. Вып. 2. 2006. С. 90-101.

Статья поступила в редакцию 15 июня 2010 г.