CC BY
10
DOI: 10.15593/2224-9982/2020.62.04 УДК 534.015:62-755
Н.Н. Зайцев1, Д.Н. Зайцев1, А.А. Макаров2, Д.А. Минеев1
1Пермский национальный исследовательский политехнический университет, Пермь, Россия 2Научно-производственное объединение «Искра», Пермь, Россия
ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ДИНАМИКИ ГИБКОГО РОТОРА С ДВУМЯ ШАРОВЫМИ АВТОБАЛАНСИРАМИ
Шаровые автобалансирующие устройства могут осуществлять компенсацию изменений неуравновешенности «на ходу» только у роторов, работающих на сверхкритических частотах вращения. Для автоматического уравновешивания таких роторов, классифицируемых как гибкие роторы, необходимо несколько автобалансиров, располагаемых в разных сечениях вала. Это обусловливает необходимость учета при исследованиях динамики ротора с автобалансирами его изгиб-ных колебаний, что особенно важно при проектировании реальных роторов. Ввиду сложности экспериментальных исследований таких роторов в статье рассматривается методика прямого численного моделирования динамики системы гибкий ротор-опоры-автобалансиры. Ее методологической основой является использование дискретной многомассовой модели ротора, эквивалентной по динамическим характеристикам реальному ротору, и уравнений динамики системы дискретный ротор-опоры-автобалансиры, полученных в прямой форме записи. Для определения дискретных масс и матрицы коэффициентов влияния жесткостей сечений ротора предполагается использование расчетов для конечно-элементной модели реального ротора в существующих программных комплексах инженерного анализа. Полученная методом Лагранжа математическая модель динамики системы учитывает нестационарность частоты вращения ротора, воздействие сил тяжести и трение качения шаров в обоймах автобалансиров. Верификация математической модели осуществлена воспроизведением опубликованных данных с использованием вычислительной модели для двухопорного однодискового трехмассового ротора с двухшаровым автобалансиром. Для четырехмассового ротора с двумя двухшаровыми автобалансирами приводятся результаты численного моделирования динамики на режимах разгона, установившегося вращения и торможения. Показано, что для рассматриваемой системы на режиме установившегося вращения имеет место лишь частичная автобалансировка, в том числе после ступенчатого увеличения дисбаланса.
Ключевые слова: шаровое автобалансирующее устройство, динамика гибкого ротора, дискретная модель ротора, двухшаровой автобалансир, результаты численного моделирования, нестационарная частота вращения, сила тяжести, трение качения.
N.N. Zaytsev1, D.N. Zaytsev1, A.A. Makarov2, D.A. Mineev1
Perm National Research Polytechnic University, Perm, Russian Federation 2Research and Production Association "Iskra", Perm, Russian Federation
NUMERICAL SIMULATION OF THE DYNAMICS OF A FLEXIBLE ROTOR WITH TWO BALL AUTO-BALANCERS
Ball auto-balancing devices can to compensate changes of unbalance "on the move" only for rotors operating at supercritical speeds. For automatic balancing of such rotors, classified as flexible rotors, several auto-balancers located in different cross sections of the shaft are necessary. This makes it necessary to account bending fluctuations on studies of dynamics of the rotor with auto-balancers, that is especially important in the design of the real rotors. In view of the complexity of experimental studies of such rotors in the article the method of direct numerical simulation of the dynamics of the flexible rotor system - supports - auto-balances is considered. The methodological basis of this method is the use of a discrete multi-mass rotor model, which is equivalent in dynamic characteristics to a real rotor, and also the equations of dynamics of the system discrete rotor - supports - auto-balancers, obtained in the direct form of recording. For definition of discrete masses and a matrix of coefficients of influence of stiffness of rotor cross-sections it is supposed to use calculations for finite-element model of a real rotor by existing software complexes of the engineering analysis. The mathematical model of the system dynamics obtained by the Lagrange method takes into account the non-stationarity of the rotor rotation speed, the influence of gravity and the rolling friction of the balls in the auto-balancer cages. Verification of the mathematical model was performed by reproducing the published data using a computational model for a two-support single-disk three-mass rotor with a two-ball auto-balancer. For a four-mass rotor with two two-ball auto-balancers, the results of numerical simulation of dynamics for the modes of acceleration, steady-state rotation and deceleration are presented. It is shown that for the system under consideration, only partial auto-balancing takes place in the steady rotation mode, including after a stepwise increase of the imbalance.
Keywords: ball auto-balancing device, flexible rotor dynamics, discrete rotor model, two-ball auto-balancer, numerical simulation results, non-stationary rotation speed, gravity, rolling friction.
Введение
Применение шаровых автобалансирующих устройств (АБУ) для компенсации «на ходу» изменений неуравновешенности роторов возможно только для роторов, имеющих сверхкритические рабочие частоты вращения [1]. Для таких роторов, классифицируемых как гибкие, характерны изгибные колебания. При этом для компенсации динамической неуравновешенности ротора требуется несколько АБУ, располагаемых в различных сечениях ротора [2]. Это обусловливает необходимость учитывать при исследованиях динамики системы ротор-опоры-АБУ формы изгибных колебаний ротора и взаимовлияние жесткостных свойств его сечений.
По исследованиям динамики роторов с несколькими шаровыми АБУ имеется много публикаций, например работы [2-9], в большинстве которых рассматривается построение математических моделей и анализ условий возникновения и устойчивости режимов автобалансировки роторов с различной неуравновешенностью. Так, в работе [3] приводится аналитическое и численное исследование динамики сплошного абсолютно жесткого ротора в упруго-демпферных опорах с двумя двухшаровыми АБУ. Авторами отмечается, что на процесс балансировки влияют как конструктивные и режимные параметры, так и начальные условия численного моделирования. Также указано, что полная балансировка в их исследованиях имела место только при скоростях вращения выше второй критической частоты.
Модели гибких роторов с несколькими шаровыми АБУ исследованы в работах [7-9]. В работе [7] для обобщенной модели ротора в виде невесомого вала с произвольным расположением конечного числа несбалансированных точечных масс и шаровых АБУ его из-гибные деформации задаются аппроксимацией на основе методов сопромата линий прогиба на участках между точечными массами и АБУ. По результатам аналитического исследования и численного моделирования авторы констатируют, что для получения частичной балансировки скорость вращения должна быть больше первой критической час-
тоты, шаровые АБУ стремятся минимизировать вибрации только в случае, когда параметры системы находятся в области устойчивости равновесного состояния, полная балансировка может достигаться при расположении несбалансированных масс в плоскостях АБУ. Работы [8, 9] посвящены математическому моделированию динамики гибких роторов с двумя и более АБУ с исследованием конструктивных и режимных условий процессов автобалансировки. В работе [8] модель гибкого ротора представлена дискретным многомассовым ротором с учетом матрицы коэффициентов же-сткостей сечений с точечными массами, в работе [9] для гибкого ротора используется уравнение упругой линии вала постоянного сечения при установившемся вращении ротора. В большинстве указанных выше работ при исследованиях не учитываются силы тяжести и невязкого сопротивления.
В данной статье рассматривается методология построения математической модели гибкого ротора с несколькими шаровыми АБУ с учетом динамических коэффициентов взаимовлияния сечений ротора, сил тяжести, а также сил вязкого и контактного трения шаров в обоймах АБУ. Для четырехмассового ротора с двумя двухшаровыми АБУ приводятся результаты численного моделирования динамики на режимах разгона, установившегося вращения и торможения.
Методологический подход к моделированию динамики гибкого ротора
В качестве методологической основы моделирования гибкого ротора выбрано представление континуального ротора дискретной многомассовой моделью [10], эквивалентной по динамическим свойствам его конечно-элементной (КЭ) модели. Для построения дискретной модели используется методика, описанная в работах [11, 12], которая предполагает для определения дискретных масс и коэффициентов влияния же-сткостей проведение расчетов в программном комплексе (ПК) АШУ8.
Далее без потери общности используется четырехмассовая дискретная модель, расчетная схема которой показана на рис. 1.
,м2
^ ^—о г2 У2 < >
г3
я, г
Рис. 1. Расчетная схема четырехмассового дискретного ротора
Согласно указанной методике для заданного расположения масс М, их величины определяются по данным из расчетов в АК8У8, исходя из требуемого совпадения учитываемых собственных частот и форм колебаний для конечно-элементной и дискретной моделей ротора. Элементы матрицы Д эквивалентных коэффициентов податливостей с учетом принципа взаимности рассчитываются по формулам, полученным в соответствии с расчетной схемой на рис. 1:
1 I - г А,, =—; Д12 = Д21 =-1;
I - г3
Д13 =Д31 =—; Д14 = Д41 =0;
с11
Д = * (£-^2)2 г22 . Д =Д =„
Д22 ®22 + „2 + л2 ' Д23 Д32 "23 +
. - г2)(^ - гз) + Д = Д =
+ с/ + с/2; Д24 Д42 с/'
(1)
д =5 +(1 - г3)2 н—г^2—.
Д33 °33 + с/ + с4 £2'
А =А = А = 1
Д34 = Д 43 = л ; А44 =
где с1, с4 - задаваемые в опорных точках коэффициенты жесткости; 522, 523 - коэффициенты податливости, определяемые в ПК АК8У8 для КЭ модели в межопорных точках в процессе вычисления масс.
Матрица г коэффициентов влияния же-сткостей находится обращением матрицы по-датливостей: г = А-1. Заметим, что коэффициенты жесткостей также могут быть непосредственно определены по КЭ модели в ПК
Построение математической модели
В основу построения математической модели положены следующие допущения:
1) континуальная модель реального ротора представлена эквивалентной дискретной четырехмассовой моделью. Соответственно, расчетная схема представляет собой (рис. 2) невесомый упругий вал, на котором располагаются четыре сосредоточенные массы: две в опорах (точки 1 и 4) и две между опорами (точки 2 и 3);
2) считается, что в точках 2 и 3 располагаются соответственно два шаровых АБУ, массы дисков которых с обоймой без шаров включены в массы этих точек. Центры дисков АБУ располагаются на оси вала.
3) с учетом малости допустимых прогибов реального ротора принимается, что при колебаниях ротора каждая масса перемещается в соответствующей плоскости, перпендикулярной оси расточки подшипников опор, рассматриваемой в качестве оси вращения идеально уравновешенного ротора. Таким образом, не учитывается влияние гироскопических моментов на вращающийся ротор. Соответственно, считается, что диски АБУ совершают плоско-поступательные движения;
4) в каждом автобалансире имеется по два шара с одинаковыми массоинерцион-ными характеристиками. Шары без проскальзывания катятся по внешней стороне как бы параллельных, то есть не сталкиваясь, кольцевых дорожек обоймы АБУ, на них действуют силы вязкого сопротивления среды, трения качения и тяжести.
5) полагается, что приводной двигатель ротора обладает неограниченной мощностью.
г
Рис. 2. Расчетная схема: Ок, к = 1...4, - центры сечений вала с координатами по оси О2 соответственно Мк, ак - сосредоточенные в сечениях Ок массы и векторы их эксцентриситетов; ук - фазовый угол к-го вектора дисбаланса Мк ак, ук = у + у0к, где у - фазовый угол вращения, у0к - начальный фазовый угол; гк, хк, ук - соответственно вектор смещения и координаты центра Ок в неподвижной СК; фк„ к = 2, 3, I = 1, 2, - угловая координата г-го шара в обойме к-го АБУ, фкг = акг + ук, где акг - угловая координата 1-го шара во вращающейся СК Окх'у'; Як - радиус орбиты движения центров шаров Вкг в обойме
Соответственно, нестационарные и стационарные режимы вращения ротора определяются задаваемыми законами изменения частоты вращения;
6) виды неуравновешенности ротора задаются величинами и угловыми координатами эксцентриситетов центра масс ротора (далее используется задание векторов дисбаланса межопорных масс);
7) крутильные колебания ротора не рассматриваются;
8) ротор и опоры считаются изотропными с не зависящими от частоты вращения коэффициентами жесткости и демпфирования.
Выбраны следующие системы координат (СК), аналогичные примененным в работах [1, 13]:
1. Неподвижная СК Оху2 с началом О на оси расточки опорных подшипников. Оси Ох и Оу в плоскости, перпендикулярной оси подшипников. Ось О2 направлена по оси подшипников, образуя правую систему координат, и является осью плоскостного прецесси-
онного движения рассматриваемых центров Ок, к = 1.4, сечений вала с координатами по оси О2 соответственно 2к.
2. Подвижные невращающиеся СК Окх'0у0 с началом Ок в центре к-го сечения вала. Плоскости Окх'0у'0 в своем движении всегда параллельны плоскости Оху неподвижной СК. Оси подвижных СК всегда сонаправлены соответствующим осям неподвижной СК.
3. Вращающиеся СК (ВСК) Окх'у' с началом Ок в центре к-го сечения вала. Оси Окх' направлены по вектору заданного в данном сечении эксцентриситета центра масс. Оси Окх' и Оку' вращаются вокруг начала координат с направлением и частотой собственного вращения ротора: га = у , где у - фазовый угол поворота ротора. Фазовые положения вращающихся СК определяются их фазовыми углами поворота ук = у + у0к, где у0к - начальный фазовый угол, отсчитываемый от положительной полуоси оси Окх'0. Соответственно,
в СК Okx'y' диски автобалансиров являются неподвижными (невращающимися).
Координаты используемых СК связаны следующими соотношениями:
- для Okx'y' и Okx'oy'o
x0 = x' cos yk - y' sin yk, y'o= x'sin yk + y'cos yk;
x' = x' cos yk + y'sin yk , y' = -x0sin yk + y'0 cos yk;
- для ox y0 и Oxy
x = xk + x0, У = Ук + У0,
Для построения математической модели динамики используются уравнения Лагранжа 2-го рода:
í дт \
dt
dq
дт ди
Ч дЧг
dW
dq,
= Q, ,
где T и U - соответственно кинетическая и потенциальная энергии системы; W - диссипа-тивная функция, учитывающая влияние сил вязкого сопротивления; q и q¿ - соответственно обобщенная координата и ее производная по времени; Qq - обобщенная сила.
В качестве вектора обобщенных координат рассматривается
q = (^ Уl, x2, У2, Х3, У3, x4, У4, Ф21, Ф22, Ф31, Фз2 ) ,
где (см. рис. 2): xk, yk, k = 1.4, - координаты в неподвижной СК центров сечений вала Ok; Ф k1, Ф k2, k = 2, 3, - угловые координаты в СК Okx0 y0 соответственно 1-го и 2-го шара в двух АБУ, расположенных в сечениях 2 и 3.
Кинетическая энергия системы с учетом энергии качения шаров и в предположении отсутствия эксцентриситетов центров масс в сечениях O1 и O4 имеет вид
T = 0,5M1 x2 + 0,5M1 y2 + 0,5M4 x42 + + 0,5Ma y^ + 0,5JCTy2 + 0,5 Jc 4 y2 + + 0,5M2(x2 - a2y sin y2)2 + + 0,5M2(y2 + a2y cos y 2)2 + 0,5 Jc2y2 +
+ 0,5M3(x3 - a3y sin y 3)2 + + 0,5M3( y3 + a2 y cos y 3)2 + + 0,5 Jc3y2 + 0,5m2 (x2 - Я2ф21 sin ф21 )2 + + 0,5m2 (y2 + R2 ф21 cos ф21)2 +
+ 0,5 J2[ R2 (Ф21 -y) - r y ]2 +
+ 0,5m2(x2 - R2ф22sin ф22)2 + + 0,5m2( y2 + R2ф22cos ф22)2 +
+ 0,5J[R2^22 -y) -Г2y]2 +
+ 0,5m3(x3 -R3ф31 sinф31)2 + + 0,5 m3( y3 + R3ф31 cos ф31)2 +
+ 0,5 J3[ R3 (ф31 -y) - Г3 y ]2 +
r3
+ 0,5m3(x3 -R3ф32 sinф32)2 + + 0,5m3(y3 + R3ф32 cos ф32)2 +
+ 0,5 J3[Rз(фз2 -y)-Г3y]2,
где (см. рис. 2): Ы\, М4 - массы в сечениях О1 и О4 (массы цапф) без эксцентриситетов; М2, М3 - массы без шаров в сечениях О2 и О3 с заданными эксцентриситетами соответственно а2 и а3; Jc4 и Jc2, ^с3 - полярные моменты инерции соответственно цапф и дисков с АБУ относительно своих центральных осей; у -угловая скорость (частота) собственного вращения ротора; у2, у3 - фазовые углы поворота дисбалансов Мкак, к = 2, 3, вокруг своих центров Ок; т2, т3 - массы шаров в АБУ, расположенных соответственно в сечениях О2 и О3 (в каждом АБУ два одинаковых шара, т.е. т21 = т22 = т2 и т31 = т32 = т3); Я2, Я3 - радиусы орбит движения центров масс (геометрических центров) шаров по обойме АБУ, расположенных соответственно в сечениях О2 и О3; J2, J3 - моменты инерции шаров в АБУ, расположенных соответственно в сечениях О2 и О3 ^21 = J22 = J2 и Jзl = Jз2 = Jз); Г2, Г3 - радиусы шаров в АБУ, расположенных соответственно в сечениях О2 и О3 (Г21 = Г22 = Г2 и Г31 = Г32 = Г3); ф2г, ф3г, г = 1, 2, - угловые координаты шаров соответственно в СК О2х'у', и О3х'у'.
Потенциальная энергия и системы определяется как
x0 = x - xk, y0 = у - yk .
3
U = Ui + U2,
где
Ui = 0,5ZZ rAXi + 0,5Hri1JtJ,
- потенциальная энергия от упругих сил системы [14, с. 82]; rki - коэффициенты влияния жесткостей системы, удовлетворяющие принципу взаимности (rki = r ik);
U2 = Mi gVi + M4 g^ + Ms 2 gJ2 +
+M2 ga2 sin y 2 + m2 gR2 sin ф21 +
+ m2 gR2 sin Ф22 + MS 3 gy3 +
+M3 ga3 sin y 3 + m3 gR3 sin ф31 + m2 gR3 sin ф32
- потенциальная энергия от сил тяжести цапф и дисков с шаровыми АБУ [13], где MS2 = M2 + 2m2, MS3 = M3 + 2m3.
Диссипативная функция Ж вязкого сопротивления в системе принимается в виде
k = 2, 3; i = 1, 2,
Ж = Wi + W2,
где
где k - номер сечения с АБУ; i - номер шара; цк - коэффициент трения качения шаров в k-м АБУ.
Подстановкой в уравнения Лагранжа выражений для T, U, W и Qqi после преобразований получается система нелинейных дифференциальных уравнений динамики че-тырехмассового ротора с двумя двухшаровы-ми АБУ, имеющая 24-й порядок:
i Xj P11 Г1 — P12 ^2 Pi 3 X3
Pl4X4 _ Г12X2 _ Г13X3 _ Г14X4;
M1 У +P11 У + Г11У —-P12y2 -P13У "
" P14У4 -Г12У2 -Г13У3 -Г14У4 -M1 g;
MS 2 X2 +P22 X2 + r22 x2 — — M2a2 (у sin у2 + у2 cos у2) + + m2 R2^21sin ф21 + ф21 cos ф21) + +m2 R2 (ф22 sin ф22 + ф22 cos ф22) -
- P12 X1 - P 23 X3 - P 24 X4 - Г12 X1 - Г23 X3 - Г24 X4 ;
к=1 ,=1 к=1 ,=1
- диссипативная функция Релея для описания вязкого сопротивления при колебаниях ротора [14, с. 132]; рк - диссипативные коэффициенты, удовлетворяющие принципу взаимности (рй = р,*);
^ =0,5]Г[6А2(фк1 -у)2 +ькяК^к2 -у)2]
к=2
- диссипативная функция, учитывающая влияние на движение шаров вязкого сопротивления среды [13]; Ьк, к = 2, 3, - коэффициенты вязкого сопротивления движению шаров в обойме АБУ.
Обобщенные силы в рассматриваемой системе принимают вид [13]
Q*t — Qyt — ^ k = L--4;
Q<tk,— -Rk—m (Rk ф«- xkcos фи -
-yk sin фkI.) sgn(фkI■ -y),
MS 2 У2 +p22 y2 + r22 У2 —
— -M2a2 (y cos у 2 - у2 sin у 2) -- m2R2 (ф21 cos ф21 - ф^ sin ф21) --m2R2 (ф22 cos ф22 - ф^2 sin ф22) -
- P12 y1 - P 23 y3 - P 24 y4 - Г12 У1 -
- Г23 y3 - Г24 у4 - MS 2 g;
(m2 + Jt)R2ф21 + Ь2R2 (ф21 -у) — Г2
— m2 (X2 sin ф21 - y2 cos ф21) -
- m2g cos ф21 + J (R2 + Г2)у -
- m2 —( R2 ф21 - X2 cos ф21 -Г2
- y2sin Ф2l)sgn(ф2l-у); (2)
(m2 + J7)R2ф22 + Ь2R2 (ф22 - у) — Г2
— m2 (x2 sin ф22 - y2 cos ф22) - m2g cos ф22 + + —2-(R2 + r2)у-m2 — (R2ф^2 -X2cosф22 -
- y2 sin ф22^п(ф22 - у);
2
k
S3x3 + P33x3 I Г33 x3 — =M3a3(y sin y3 + y2 cos y3) + + m3 Д,(ф31 sin ф31 +ф21 cos ф31) +
+ m3R3 (ф32 sin Ф32 + ф32 cos ф32 ) -
- P13 x1 - P23x2 - P34x4 - r13x1 - Г23x2 - Г34 x4 ;
MS 3 y3 +p33 y3 + r33 У3 =
= -M3a3 ( y cos y3 - y2 sin y 3) -
- m3 R3 (ф31 cos ф31 - ф^ sin ф31) -
- m3R3 (ф32 cos ф32 - ф^2 sin ф32) -
- р13y1 -P23y2 -p34y4 - r13У1 -- Г23У2 - r34У4 - MS 3 g;
(m + J2-) RзФзl + ¿3 Rз(Фзl -y) =
r3
= m3 (x3 sin ф31 - y3 cos ф31) - m3 g cos ф31 + +—3(R3 + r3) y-m3 —(R,^ - x3cos ф31 -
r3 r3
- У3 sin Фзl)sgn(фзl -y);
(m3 + 4-)RзФз2 + ¿3Rз(Фз2 -y) =
Г3
= m3 (x3 sin ф32 - y3 cos ф32) - m3g cos ф32 +
+JF(R3 + r3)y - m3 —(Rзф22 - x3cosФ32 -
r3 r3
- y3 sin Ф32 ) sgn(ф32 -y);
M- 4x4 + P44x4 + r44x4 = -P14xl - P24x2 -- P34 x3 - Г14 xl - Г24 x2 - Г34 x3;
M4У4 + P44У4 + Г44У4 = -P14У1 - P24У2 -- P34У3 - Г14У1 - Г24У2 - Г34У3 -M4g.
Здесь (см. рис. 2)
Фи =ak, +y k, k=2, 3; ■ = 1, 2 (3)
где ak - угловая координата i-го шара во вращающейся СК Okx'y';
y(t) = jra(x)d х
- фазовый угол вращения ротора (см. описание вращающейся СК); ю - скорость собственного вращения ротора, задаваемая желаемым законом изменения, ш = у [13].
Ниже в расчетах при установившемся режиме задавалась постоянная скорость ш = const, а режимы разгона и торможения производились с постоянным ускорением ш = const.
Верификация вычислительной модели
Вычислительная модель динамики гибкого ротора с АБУ получается приведением методом исключений системы уравнений (2) к форме Коши с учетом уравнения (3). При этом уравнения динамики для межопорных масс с АБУ (см. M2 и M3 на рис. 1) имеют такой же вид, как в работе [13]. Численное интегрирование системы в форме Коши осуществляется методом Рунге - Кутты 4-го порядка с использованием описанной в работе [13] методики.
Верификация вычислительной модели осуществлялась воспроизведением результатов численного моделирования одномассового симметричного ротора с АБУ (прототип) из работ [15] и [13] с помощью математической модели динамики для трехмассового ротора (две массы Mi и M3 в опорах и одна M2 с АБУ посередине пролета). Матрица эквивалентных податливостей в этом случае рассчитывалась соответственно формулам (1) при z2 = 0,5£ . Массы задавались следующим образом: M2 = 0,1 кг, что соответствовало прототипу, а M1 = M3 = 0,000 05 кг. Моделировались два случая: первый при жестких опорах с с1 = с3 = = 500 000 Н/м, второй при податливых опорах с с1 = c3 = 500 Н/м. В обоих случаях податливость вала в сечении с АБУ принималась как у прототипа 522 = 0,001 м/Н. Для определения критических частот вращения рассматриваемого трехмассового ротора использовалось частотное уравнение из работы [10] (см. уравнение (2.20)). Соответственно, для первого случая они имели значения pl = 99,95, p2 = 100 000, p3 = 100 050, для второго - pl = 70,7, p2 = 3162,3, p3 = 4472,7. У прототипа единственная критическая частота pl = 100 рад/с. При моделировании конструктивные и режимные параметры соответствовали прототипу (см. работу [13]), за исключением того, что в первом случае для массы M2 задавались эксцентриситет и коэффициент демпфирования, как у прототипа, во втором - этот эксцентри-
ситет принимался 150 мкм, а коэффициенты демпфирования задавались только в опорах.
Результаты численного моделирования приведены на рис. 3-5. На рис. 3 и 4 показано изменение модуля вектора прогиба (смещения)
центра соответствующего сечения, гк = ^хк + у2к ,
к = 1, 2, 3. На рис. 5 показано изменение угловых координат (3) шаров во вращающейся СК.
Графики на рис. 3, а; 4, а и 5 в целом идентичны аналогичным графикам прототипа [13, 15]. Графики прогибов на рис. 3 и 4 соответствуют теории роторной динамики [13]. Как следует из графиков на рис. 3, б и рис. 4, б, вычислительная модель адекватно реагирует на изменение жесткостных и демпфирующих свойств в опорах.
Рис. 3. Изменение модуля вектора прогиба при жестких опорах: а - в сечении 2, б - в сечениях 1 и 3; ю - частота вращения
Рис. 4. Изменение модуля вектора прогиба при податливых опорах: а - в сечении 2, б - в сечениях 1 и 3; ю - частота вращения
Рис. 5. Изменение угловых координат шаров во вращающейся СК: а - при жестких опорах; б - при податливых опорах
Таким образом, проведенные для трех-массового ротора с двухшаровым АБУ расчеты позволяют считать разработанную вычислительную модель динамики дискретного ротора приемлемой для проведения исследований численным моделированием динамики гибкого ротора с шаровыми АБУ.
Численное моделирование динамики четырехмассового ротора с двумя двухшаровыми АБУ
Моделирование проводилось для горизонтального модельного ротора, эквивалентная дискретная четырехмассовая модель которого приведена в работе [16]. Соответственно, использовались следующие исходные данные для вычислительной модели (см. рис. 1 и 2): M1 = M4 = 22 кг, M2 = 34 кг, M3 = 77 кг, z1 = 0, z2 = 0,4595 м, z3 = 1,1945 м, z4 = l = 1,822 м,
522 = 0,778-10-8 м/Н, 5зз = 0,1091-10-7 м/Н,
523 = 5з2 = 0,7679-10-8 м/Н, cx = С4 = 108 Н/м. Для двух одинаковых двухшаровых АБУ, располагаемых в сечениях с массами M2 и M3, задавались следующие параметры: R2 = R3 = 0,08 м, m2 = m3 = 15,9 г, r2 = r3 = 7,87 мм, J2 = J3 = = 3,95-10-7 кгм2, b2 = b3 = 0,05 Нс/м, Ц2 = M3 = = 0,000 05 м. В этих же сечениях задавалась неуравновешенность параметрами a2 = a3 = = 40 мкм и у02 = у03 = 0. Начальные значения переменных, кроме положения шаров, задавались нулевыми. Начальные положения шаров в соответствующих вращающихся СК Okx'y' принимались равными a2i = a3i = -93o и a22 = a32 = -87o. Гравитационное ускорение g = 9,81 м/с2.
Элементы матрицы коэффициентов влияния жесткостей, найденные обращением матрицы эквивалентных податливостей, рассчитанной по формулам (1), имели значения (Н/м): Г11 = 2,18-108, Г12 = Г21 = -2,13-108, Г13 = Г31 = = 1,18-108, Г14 = Г41 = -2,38-107, Г22 = 4,21 -108, r23 = Г32 = -2,96-108, Г24 = Г42 = 8,81-107, Г33 = 3,00-108, Г34 = Г43 = -1,22-108, Г44 = 1,58-108. Демпфирование задавалось только в опорах Р11 = р44 = 3,5-104 Нс/м.
Расчеты производились методом Рунге -Кутты 4-го порядка в программной среде Processing 3 с возможностью построения диа-
грамм и анимации. Численное интегрирование выполнялось с шагом 0,000 01 с (частота дискретизации 100 000 Гц) при выводе на печать с шагом 0,0006 с (частота дискретизации 1666,667 Гц). Частота установившегося вращения продолжительностью 48 с составляла 162,5 Гц. Продолжительность разгона и торможения принималась одинаковой, равной 12 с. На режиме установившегося вращения в момент t = 36 с увеличивался эксцентриситет a3 в 1,3 раза. Вычисления производились для ротора без АБУ (без шаров) и с АБУ (с шарами). Результаты расчетов представлены на рис. 6-9.
На рис. 6 и 7 показаны графики изменения модулей векторов прогиба (смещения)
rk = Vхк + Ук , k = 1.. .4, центров сечений Ok на режимах разгона, установившегося вращения и торможения. Согласно графикам АБУ не обеспечивают полной автобалансировки на установившемся режиме, но значительно снижают амплитуду колебаний ротора, в том числе после увеличения одного из дисбалансов. Можно также видеть, что наличие АБУ при прохождении критической частоты во время разгона приводит к некоторому увеличению максимума амплитуды колебаний, а при торможении, наоборот, к его уменьшению.
Для оценки влияния автобалансиров на динамическое поведение рассматриваемого ротора с использованием кинетостатики Даламбе-ра и координатных преобразований рассчитывались проекции главного вектора (P) и главного момента (M) сил инерции относительно мгновенного центра масс ротора. На рис. 8 и 9 показаны изменения проекций на плоскость Oxy соответственно P и M, определяемых как
Pxy =VPx2 + Py и Mxy =VMTM. Вычисления производились с учетом и без учета сил тяжести, влияние которых для данного ротора незначительно.
На рис. 8 и 9 видно, что на установившемся режиме вращения величины главного вектора и главного момента сил инерции у ротора с автобалансирами меньше, чем при их отсутствии, что свидетельствует об уменьшении автобалансирами статической и момент-ной неуравновешенностей ротора.
Рис. 6. Смещение центров сечений 02 и 03 для ротора: а - без АБУ; б - с АБУ
20200828200150
Модуль прогиба, частота вращения
2 1,00 я" 0,80
тт,-и
^ 0,60
I 0,40
о
0,00
~Т ч- \
t~ г
-i —h Т—
Р Л- h te-_
0,20 0,16
20200828185700
Модуль прогиба, частота вращения 0 1,00""""""
0,43
г-,2
0,00 0,14 0,29
U с, Ю'
in3d019_LeykihYeBaL]_ga l2_40mkrNLT_180826-.xls
а
0,58 0,72
0,12 |
0,08 | О
0,04 Я
О
0,00 У
У
0,80
0,60
I 0,40
I 0,20 g
I 0,00
/ /
/ \
/
J
!
/
/ 1 \
и 1 V
0,00 0,14 0,29
0,43
0,20 -0,16 £
0,12 | a
0,08 g, и
0,04 Я
О
0,00 % 0,58 0,72 ^
t, с, 10"
in3d019 LeyMYeBall g al2 40mkrNU'J8082<j'.sls б
Рис. 7. Смещение центров сечений 01 и 04 для ротора: а - без АБУ; б - с АБУ
Рис. 8. Изменение проекции на плоскость Оху главного вектора сил инерции: а - без АБУ; б - с АБУ, при расчетах с учетом (О) и без учета (N00) сил тяжести
20200828200150
20200828185700
0,20 ■о 0,16
о '
Я 0,12 К
^ 0,08 ^ 0,04 0,00
Модуль (от X, У-проекций) главного вектора сил (Оху), частота вращения
/ \
/ {
\
/
/ ~т {V
т г \
о о
0,20 2 я
0,16
3
0,12 я
' и
0,08 | ' Рч
ш
0,04 Й о
н
0,00 0,14 0,29 0,43 0,58 0,72 1, с, 102
¡пЗй019_ЬеукШУеВа1^ а12_40ткгШ _180826-.xls
а
0,00 Й
0,20 о 0,16
о '
а 0,12
К
_1 0,08 ^ 0,04 0,00
Модуль (от X, У-проекций) главного вектора сил (Оху), частота вращения
т
{ \
ч- л—
п Л
0,20 2 0,16 ^ 0Л2 |
0,08 I
га
0,04 £
0,00 0,14 0,29 0,43 0,58 0,72 1, с, Ю2
¡пЗй019_ЬеукШУеВа1^ а12_40ткгШ 180826 .хЬ
0,00
Рис. 9. Изменение проекции главного момента на плоскость Оху для ротора: а - без АБУ; б - с АБУ, при расчетах с учетом (О) и без учета (N00) сил тяжести
б
Рис. 10. Скриншоты анимации: а - момент максимальной амплитуды при разгоне ($ = 9 с); б - режим установившегося вращения при t = 59 с
а
На рис. 10 приведены скриншоты анимации в моменты вращения ротора с АБУ на режимах разгона и установившегося вращения. Здесь для сечений О2 и О3 разноцветными стрелками обозначены направления радиус-векторов: черные - для заданных эксцентриситетов масс (осей Окх'), красные - для прогибов, зеленые - для мгновенных центров масс сечений с шарами, синие - для мгновенных дисбалансов, создаваемых только шарами. Для перечисленных в столбцах параметров указаны текущие значения их угловых координат (в градусах в соответствующих СК) и амплитуды (в микрометрах). В правом верхнем углу располагается диаграмма текущих значений координат смещений центров сече-
ний и соответствующего им модуля прогиба (зеленая ломаная линия).
На рис. 10, а можно видеть, что на режиме разгона при прохождении критической частоты вращения шары находятся на «тяжелой» стороне дисков с автобалансирами, т.е. по одну сторону с эксцентриситетами масс данных сечений, увеличивая тем самым амплитуду колебаний ротора. На установившемся режиме вращения (см. рис. 10, б) шары находятся уже на «легкой» стороне, т.е. диаметрально относительно эксцентриситетов масс данных сечений, что приводит к уменьшению амплитуды колебаний [1]. Следует отметить, что в процессе расчета на установившемся режиме вращения расположение шаров относительно друг друга в обоих АБУ,
показанное на рис. 10, б, остается неизменным, хотя после скачкообразного увеличения дисбаланса угловые координаты (ак) шаров в СК 0кх'у' возрастают по абсолютной величине на несколько градусов.
Заключение
Проведенные исследования позволяют отметить следующее:
1. Поскольку шаровой автобалансир может осуществлять автобалансировку ротора только на сверхкритических частотах вращения, а для компенсации его динамической неуравновешенности требуется не менее двух АБУ, математическая модель динамики системы ротор-опоры-автобалансиры должна учитывать изгибные колебания ротора. Для этого предложено использовать дискретную многомассовую модель реального ротора, эквивалентную по динамическим свойствам его конечно-элементной модели, с определением дискретных масс и матрицы коэффициентов
влияния жесткостей на основе расчетов в программном комплексе ANSYS.
2. Построением с помощью метода Ла-гранжа математической модели динамики че-тырехмассового ротора с двумя двухшаровы-ми АБУ показана возможность получения системы обыкновенных дифференциальных уравнений в прямой форме для численного моделирования динамики системы гибкий ро-тор-опоры-шаровые автобалансиры.
3. Численным моделированием динамики системы четырехмассовый ротор-опоры-два двухшаровых АБУ установлено, что, несмотря на отсутствие полной автобалансировки ротора на установившемся режиме вращения, имеет место существенная компенсация его статической и моментной неуравновешен-ностей. Показано, что наличие АБУ обеспечивает автобалансировку данной роторной системы после скачкообразного увеличения дисбаланса, а также уменьшает максимум амплитуды при прохождении критической частоты на режиме торможения.
Библиографический список
1. Зайцев Н.Н., Зайцев Д.Н., Макаров А.А. Инженерный анализ установившихся режимов одно-дискового ротора с многорядным шаровым автобалансирующим устройством // Вестник Пермского национального исследовательского политехнического университета. Аэрокосмическая техника. - 2017. -№ 48. - C. 43-59. DOI: 10.15593/2224-9982/2017.48.05
2. Гончаров В.В., Филимонихин Г.Б. Вид и структура дифференциальных уравнений движения и процесса уравновешивания роторной машины с автобалансирами // Известия Том. политехн. ун-та. -
2015. - Т. 326, № 12. - С. 20-30.
3. Automatic two-plane balancing for rigid rotors / D.J. Rodrigues, A.R. Champneys, M.I. Friswell, R.E. Wilson // Int. J. of Non-Linear Mech. - 2008. - Iss. 43. - P. 527-541.
4. Bolton J.N. Single- and dual-plane automatic balancing of an elastically mounted cylindrical rotor with considerations of coulomb friction and gravity: Dis for the degree of Doctor of Philosophy in Engineering Mechanics. - 2010. - URL: https://vtechworks.lib.vt.edu/handle/10919/29946?show=full (accessed 19 August 2020).
5. Ковачев А.С. Балансировка динамически неуравновешенного ротора с учетом неидеальности автобалансировочных устройств // Вестник СПбГУ. - 2015. - Сер. 1, т. 2 (60), вып. 4. - С. 606-616.
6. Об ограничении точности балансировки шаровыми автобалансирами крыльчатки осевого вентилятора / Л.С. Олийниченко, В. А. Грубань, М.В. Личук, В.В. Пирогов // Восточно-Европейский журнал передовых технологий. - 2018. - Т. 1, вып. 1 (91). - С. 27-35. DOI: 10.15587/1729-4061.2018.123025
7. Ehyaei J., Moghaddam M.M. Dynamic response and stability analysis of an unbalanced flexible rotating shaft equipped with N automatic ball-balancers // J. of Sound. and Vibration. - 2009. - Vol. 321, iss. 3-5. - P. 554-571. DOI: 10.1016/j.jsv.2008.10.019
8. Research of stability and transition processes of the flexible double-support rotor with auto-balancers near support / V. Goncharov, A. Nevdakha, Yu. Nevdakha, V. Gutsul // East-Europ. J. of Enterprise Technol. -
2016. - Vol. 6, iss. 7 (84). - P. 22-27. DOI: 10.15587/1729-4061.2016.85461
9. Methods of balancing of an axisymmetric flexible rotor by passive auto-balancers / G. Filimonikhin, I. Filimonikhina, K. Dumenko, V. Pirogov // East-Europ. J. of Enterprise Technol. - 2017. - Vol. 3, iss. 7 (87). - P. 22-27. DOI: 10.15587/1729-4061.2017.101832
10. Хронин Д.В. Колебания в двигателях летательных аппаратов. - М.: Машиностроение, 1980. -
296 с.
11. Гадяка В.Г., Лейких Д.В., Симоновский В.И. Математическая модель ротора турбокомпрессора для исследования несинхронных составляющих вибрации // Компрессорное и энергетическое машиностроение. - 2010. - № 2 (20). - С. 48-50.
12. Симоновский В.И. Оценивание коэффициентов математических моделей колебательных систем: учеб. пособие. - Саарбрюккен, Германия: LAP Lambert Academic Publishing, 2015. - 100 с.
13. Зайцев Н.Н., Зайцев Д.Н., Минеев Д.А. Моделирование динамики однодискового ротора с шаровым автобалансиром на переходных и установившихся режимах вращения // Вестник Пермского национального исследовательского политехнического университета. Аэрокосмическая техника. - 2019. -№ 57. - C. 148-161. DOI: 10.15593/2224-9982/2019.57.12
14. Бидерман В. Л. Теория механических колебаний. - М.: Высшая школа, 1980. - 408 с.
15. Yoshida S., Naka T. Reduction method of residual balancing error on auto-balancer mechanism // SICE J. of Control, Measurement and Syst. Integ. - 2014. - Vol. 7, no. 3. - P. 141-146.
16. Лейких Д.В. Идентификация причин возбуждения несинхронных колебаний роторов турбокомпрессоров и способы снижения их амплитуд: дис. ... канд. техн. наук / Сум. гос. ун-т. - Сумы, 2011. - 154 с.
References
1. Zaytsev N.N., Zaytsev D.N., Makarov A.A. Inzhenernyy analiz ustanovivshikhsya rezhimov odnodiskovogo rotora s mnogoryadnym sharovym avtobalansiruyushchim ustroystvom [Engineering analysis of steady-state regimes of the single-disk rotor with multi-row automatic ball balancing device]. PNRPU Aerospace Engineering Bulletin, 2017, no. 48, pp. 43-59. DOI: 10.15593/2224-9982/2017.48.05
2. Goncharov V.V., Filimonikhin G. B. Vid i struktura differentsialnykh uravneniy dvizheniya i protsessa uravnoveshivaniya rotornoy mashiny s avtobalansirami [Form and structure of differential equations of motion and process of auto-balancing in the rotor machine with auto-balancers]. Bulletin of the Tomsk Polytechnic University, 2015, vol. 326, no. 12, pp. 20-30.
3. Rodrigues D.J., Champneys A.R., Friswell M.I., Wilson R. E. Automatic two-plane balancing for rigid rotors. International Journal of Non-Linear Mechanics, 2008, iss. 43, pp. 527-541.
4. Bolton J.N. Single- and dual-plane automatic balancing of an elastically mounted cylindrical rotor with considerations of coulomb friction and gravity: dissertation for the degree of Doctor of Philosophy in Engineering Mechanics. 2010. URL: https://vtechworks.lib.vt.edu/handle/10919/29946?show=full (Treatment date: 19/08/2020)
5. Kovachev A.S. Balansirovka dinamicheski neuravnoveshennogo rotora s uchetom neidealnosti avtobalansirovochnykh ustroystv [Balancing a dynamically unbalanced rotor taking into account the non-ideal of auto-balancing devices]. Vestnik of Saint Petersburg University, 2015, ser. 1, vol. 2 (60), iss. 4, pp. 606-616.
6. Oliynichenko L.S., Gruban V.A., Lichuk M.V., Pirogov V.V. Ob ogranichenii tochnosti balansirovki sharovymi avtobalansirami krylchatki osevogo ventilyatora [About limiting the accuracy of balancing by ball auto balancers of the axial fan impeller]. Eastern-European Journal of Enterprise Technologies, 2018, vol. 1, iss. 1 (91), pp. 27-35. DOI: 10.15587/1729-4061.2018.123025
7. Ehyaei J., Moghaddam Majid M. Dynamic response and stability analysis of an unbalanced flexible rotating shaft equipped with N automatic ball-balancers. Journal of Sound and Vibration, 2009, vol. 321, iss. 35, pp. 554-571. DOI: 10.1016/j.jsv.2008.10.019.
8. Goncharov V., Nevdakha A., Nevdakha Yu., Gutsul V. Research of stability and transition processes of the flexible double-support rotor with auto-balancers near support. Eastern-European Journal of Enterprise Technologies, 2016, vol. 6, iss. 7 (84), pp. 22-27. DOI: 10.15587/1729-4061.2016.85461
9. Filimonikhin G., Filimonikhina I., Dumenko K., Pirogov V. Methods of balancing of an axisymmetric flexible rotor by passive auto-balancers. Eastern-European Journal of Enterprise Technologies, 2017, vol. 3, iss. 7 (87), pp. 22-27. DOI: 10.15587/1729-4061.2017.101832
10. Khronin D. V. Kolebaniya v dvigatelyakh letatelnykh apparatov [Fluctuations in flight vehicle engines]. Moscow, Mashinostroenie, 1980, 296 p.
11. Gadyaka V.G., Leykikh D.V., Simonovskiy V.I. Matematicheskaya model rotora turbokompressora dlya issledovaniya nesinkhronnykh sostavlyayushchikh vibratsii [Mathematical model of the turbocharger rotor
for investigation of nonsynchronous components of vibration]. Kompressornoe i energeticheskoe mashinostroenie, 2010, no. 2(20), pp. 48-50.
12. Simonovskiy V.I. Otsenivanie koeffitsientov matematicheskikh modeley kolebatelnykh system [Evaluation of the coefficients of mathematical models of oscillating systems]. Saarbrücken, Germany: LAP Lambert Academic Publishing, 2015, 100 p.
13. Zaytsev N.N., Zaytsev D.N., Mineev D.A. Modelirovanie dinamiki odnodiskovogo rotora s sharovym avtobalansirom na perekhodnykh i ustanovivshikhsya rezhimakh vrashcheniya [Simulation of single-disc rotor dynamics with a ball autobalancer on transient and steady-state modes of rotation], PNRPU Aerospace Engineering Bulletin, 2019, no. 57, pp. 148-161. DOI: 10.15593/2224-9982/2019.57.12
14. Biderman V.L. Teoriya mekhanicheskikh kolebaniy [Theory of mechanical fluctuations]. Moscow: Vysshaya shkola, 1980, 408 p.
15. Yoshida Shuichi, Naka Teruyuki. Reduction method of residual balancing error on auto-balancer mechanism. SICE Journal of Control, Measurement and System Integration, 2014, Vol. 7, no. 3, pp. 141-146.
16. Leykikh D.V. Identifikatsiya prichin vozbuzhdeniya nesinkhronnykh kolebaniy rotorov turbokompressorov i sposoby snizheniya ikh amplitud [Identification of causes of excitation of non-synchronous vibrations of turbocharger rotors and ways to reduce their amplitudes]: Dissertation of the candidate of technical sciences, Sumy state university, Sumy, 2011. - 154 p.
Об авторах
Зайцев Николай Николаевич (Пермь, Россия) - доктор технических наук, профессор кафедры «Ракетно-космическая техника и энергетические системы» ФГБОУ ВО ПНИПУ (614990, г. Пермь, Комсомольский пр., д. 29; e-mail: znn@perm.ru).
Зайцев Денис Николаевич (Пермь, Россия) - ведущий инженер кафедры «Ракетно-космическая техника и энергетические системы» ФГБОУ ВО ПНИПУ (614990, г. Пермь, Комсомольский пр., д. 29; e-mail: rkt@pstu.ru).
Макаров Андрей Анатольевич (Пермь, Россия) - ведущий конструктор ПАО НПО «Искра» (614038, г. Пермь, ул. Академика Веденеева, д. 28, e-mail: makarovandrej@mail.ru).
Минеев Дмитрий Андреевич (Пермь, Россия) - аспирант кафедры «Ракетно-космическая техника и энергетические системы» ФГБОУ ВО ПНИПУ (614990, Пермь, Комсомольский пр., д. 29, e-mail: mda886@mail.ru).
About the authors
Nikolay N. Zaytsev (Perm, Russian Federation) - Doctor of Technical Sciences, Professor of Rocket and Space Engineering and Power Generating Systems Department, Perm National Research Polytechnic University (29, Komsomolsky av., Perm, 614990, Russian Federation; e-mail: znn@perm.ru).
Denis N. Zaytsev (Perm, Russian Federation) - Lead Engineer of Rocket and Space Engineering and Power Generating Systems Department, Perm National Research Polytechnic University (29, Komsomolsky av., Perm, 614990, Russian Federation; e-mail: rkt@pstu.ru).
Andrey A. Makarov (Perm, Russian Federation) - Lead Designer, PJSC "Research and Production Association "Iskra" (28, Academica Vedeneeva st., Perm, 614038, Russian Federation; e-mail: makarovandrej@mail.ru).
Dmitriy A. Mineev (Perm, Russian Federation) - PhD Student of Rocket and Space Engineering and Power Generating Systems Department, Perm National Research Polytechnic University (29, Komsomolsky av., Perm, 614990, Russian Federation; e-mail: mda886@mail.ru).
Получено 09.09.2020
