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

УДК 62-133.33

Д. С. Вахлярский, А. М. Гуськов

ЧИСЛЕННЫЙ АНАЛИЗ ДИНАМИКИ РОТОРА ЦЕНТРОБЕЖНОГО НАСОСА

В работе построена математическая модель ротора насосного агрегата НМ3600-230 с учетом нелинейности упругой реакции подшипников качения, нелинейности реакции щелевого уплотнения и податливости корпусов подшипников. Методом Ньютона найдены периодические движения ротора — прямая синхронная прецессия. Исследовано влияние технологических погрешностей и работы на неоптимальном режиме на характеристики прецессионных движений.

E-mail: gouskov_am@mail.ru;_vakhlyarskyDS@gmail.com

Ключевые слова: ротор насоса, щелевое уплотнение, периодическое

решение, метод Ньютона.

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

Элементы насоса, наиболее существенно влияющие на динамику его ротора, показаны на рис. 1. На расчетной схеме (рис. 2) отражено поведение элементов, которые приведены на рис. 1.

Рис. 1. Узлы насоса, существенно влияющие на динамику его ротора:

1 — муфта; 2 — вал; 3 и 9 — подшипники качения; 4 и 8 — корпуса подшипников; 5 и 7 — кольцо уплотняющее; 6 — рабочее колесо

В модели ротора вал представлен как стержень кусочно-постоянного поперечного сечения, который разбит на конечные элементы (КЭ), основанные на модели балки Бернулли — Эйлера [3].

Рис. 2. Расчетная схема ротора насоса НМ3600-230:

1 — вал, разбитый на конечные элементы; 2 — рабочее колесо; 3 — нелинейные пружины, моделирующие подшипники качения; 4 — упругие элементы, моделирующие жесткость корпусов подшипников; 5 — присоединенные массы корпусов подшипников; 6 — нелинейные упругие элементы, моделирующие щелевое уплотнение; 7 — муфта; 8 — упругие элементы, моделирующие взаимодействие муфты с ротором электродвигателя

Обобщенными перемещениями в узлах КЭ являются перемещения и углы поворота, соответствующие неподвижным осям X и У (см. рис. 2). В уравнениях движения вал представлен матрицами жесткости К5 и масс М 5 .

Рабочее колесо рассматриваем как абсолютно твердое тело массой с главными экваториальным и осевым JAW моментами инерции. В уравнения движения рабочее колесо входит посредством гироскопической матрицы Н^ и матрицы масс MW .

Аналогично рабочему колесу рассматриваем муфту как абсолютно твердое тело массой тм с главными экваториальным JEM и осевым JAM моментами инерции. Связь муфты с валом электродвигателя смоделирована линейными пружинами с матрицей жесткости Км. Исследование жесткостных параметров упругих компенсирующих муфт и их влияния на динамику рассматриваемого ротора выполнено в работах [4, 5]. Гироскопическая матрица Нм и матрица масс Мм муфты аналогичны соответствующим матрицам рабочего колеса.

Корпуса подшипников качения (КПК) в расчетной схеме (см. рис. 2) смоделированы линейными пружинами с присоединенными массами. В уравнения движения КПК входят посредством матрицы жесткости Кн и матрицы масс М н. Схема получения матриц КПК показана на рис. 3.

Рис. 3. Схема получения эквивалентного элемента корпуса подшипника:

1 — твердотельная модель КПК; 2 — КЭ-модель КПК; 3 — вычисление матрицы жесткости; 4 — вычисление первой собственной частоты; 5 — вычисление матрицы масс; 6 — эквивалентный элемент КПК

В результате реализации данной схемы получаем

1

K h

kHxx 0

_ 0 kHyy

Мн К н, (1)

где kHxx, kHyy — жесткости КПК в направлениях осей X и У соответственно; = 802 Гц — первая собственная частота парциальных

поперечных колебаний КПК.

Рабочая частота вращения ротора равна 50 Гц, парциальная частота достигает значений нескольких сотен герц. Поэтому для моделирования подшипников качения используются их статические характеристики [6]. Ротор рассматриваемого насоса установлен на роликовые сферические подшипники С 2220 К и 22220 КЕ, которые препятствуют только радиальному перемещению цапф вала, т. е. соответствуют шарнирной схеме закрепления. Такие подшипники качения (ПК) дают только радиальную реакцию. Статическая характеристика ПК получена методом конечных элементов (МКЭ) с учетом деформации роликов и колец в процессе их контактного взаимодействия. Для аппроксимации характеристик рассматриваемых ПК используется уравнение однорядного радиального шарикоподшипника [6]:

Яв = Свп¥, (2)

где Яв — модуль радиальной реакции в подшипнике; Св — коэффициент, зависящий от типоразмера подшипника; ив — модуль радиального перемещения в подшипнике.

Коэффициенты Свь и Ст для левого и правого подшипников

соответственно определяют методом наименьших квадратов. Статическая характеристика ПК 22220 КЕ показана на рис. 4.

Рис. 4. Статическая характеристика ПК 22220 КЕ

Щелевые уплотнения (ЩУ) проточной части предназначены для уменьшения перетоков жидкости через зазоры между ротором и статором под действием перепадов давления между смежными полостями многоступенчатых насосов [2].

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

Рис. 5. Перемещения в щелевых уплотнениях при смещении и повороте рабочего колеса:

1 — уплотнительное кольцо; 2 — рабочее колесо; 3 — гипотетическое положение рабочего колеса

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

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

R=Рл+РЯ, (3)

где Ри — проекция реакции на радиальное направление; р — проекция реакции на тангенциальное направление; eu — орт радиального направления; ev — орт тангенциального направления.

Рис. 6. Распределение сил в щелевом уплотнении

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

скорости д центра колеса [4, 7] вычисляем по формулам

К =-а05-ЬсК0 (5)5;

К = ЬсЬс (5,с)5;

Кс (5) = (1 + 252)/(1 -52)2,5; (4)

Ьс = (0,5^-6?) -52 )1,5,

где ас, Ьс — коэффициенты, зависящие от параметров щелевого уплотнения; 5 — радиальное перемещение центра симметрии уплотняемой поверхности колеса; 50 — номинальный зазор в уплотнении (5 = 5/50); 5 — радиальная скорость центра симметрии уплотняемой поверхности колеса; со — скорость собственного вращения ротора; в — скорость изменения угла 6 (см. рис. 6). Работа щелевого уплотнения подобна работе короткого подшипника скольжения. Исследование динамики ротора под действием сил, которые определяются выражениями, аналогичными (4), изложено, например, в работе [1].

Коэффициенты аСл, ЬСл вычисляем по формулам

а = АрлРсЬсП. Ь = \2ЪпРсьЬс^ (5)

ас = л(л , ч2 * ; Ьс = ГГТз , (5)

4 (1 + п) 50 2450

где п = 7550/Ьс — обратный коэффициент гидравлического сопротивления кольцевого дросселя [8]; Ьс — длина ЩУ вдоль образующей; Бс — радиус колеса в месте уплотнения; Ар — разность между давлением в области выхода и входа в колесо; ¡и — динамическая вязкость жидкости.

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

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

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

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

Неконцентричность цилиндрических поверхностей в ЩУ может возникнуть вследствие различных причин. Не вдаваясь в подробно-

сти появления технологических погрешностей, расцентровку в щелевых уплотнениях можно привести к некоторому начальному смещению г,0 центра симметрии рабочего колеса по отношению к центру

вращения и начальному повороту на малые углы ру и р0 вокруг поперечных осей X и Y соответственно. При этом в исходной конфигу-

^y ^У ТТТЛГ

рации возникают расцентровки eSR и eSL в правом и левом ЩУ соответственно.

При наличии технологических погрешностей вычисляем векторы полной расцентровки 5L и 6R в левом и правом ЩУ как сумму рас-центровок, связанных с перемещением rW и поворотом ф, рабочего

колеса, с начальными векторами расцентровок eyR и eyL, повернутыми в актуальную конфигурацию с помощью тензора поворота L1 (cotk). Векторы имеют вид

5L = rW - L1G<Pw Х k + L1 (Ct k )'(( - L1GфW X k);

. . = J ! (6)

5 R = ?W + LlGфw X k + Lj (cot k ) • (W + LlGфW X k),

где k — орт оси Z; L1G — расстояние от центра симметрии колеса до середины уплотнения.

Выражение для тензора поворота L1 (cotk) имеет следующий вид [9]:

L1 (cot k) = cos (cot)E + (1- cos (cot))kk + sin (cot)EXk. (7)

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

Поскольку щелевых уплотнений два, то их воздействие на колесо

приводим к суммарной силе Rg и моменту Mg , которые вычисляем

по формулам

R G = rGL+ rGR ; (8)

MG = L2gkx(rGR - RGL),

где RgL , RgR — реакция правого и левого уплотнений соответственно; L2g — расстояние от центра симметрии колеса до проекции точки приложения реакции правого уплотнения на ось Z (см. рис. 5).

В данной работе исследуем динамику горизонтального ротора центробежного насоса при действии на него силы тяжести тW§,

инерционной нагрузки ^ и гидродинамической силы со стороны

перекачиваемой жидкости ¥ь, в качестве которой рассматриваем

воду. На рис. 7 показано рабочее колесо с действующими на него нагрузками.

Рис. 7. Нагрузки, действующие на рабочее колесо

Для пяти различных центральных угловых расположений рабочего колеса в потоке жидкости выполнены гидравлические расчеты и

найдены проекции ¥ы и Еку гидродинамической силы ¥к на оси X и

У [10—12]. Интерполяция зависимости этих проекций от времени показана на рис. 8. Поскольку колесо имеет семь лопастей, то период гидродинамической силы равен Тн = 2п/ 7с.

Рис. 8. Интерполяция проекций гидродинамической силы

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

Уравнение движения рассматриваемой расчетной схемы имеет

вид

Мх + (Б + Н)х + Кх = Г () + Я (х, х), (9)

т

где х = Ц,..., ик, ук ,фхк ,фук,..., фп } — вектор узловых перемещений;

М — матрица масс модели; К — матрица жесткости модели; Н — матрица гироскопических сил; Б = аК — матрица сил диссипации; Г (() — вектор узловых сил модели; Я (х, х) — вектор нелинейных

реакций подшипников и щелевых уплотнений. Матрицу диссипации Б = аК принимаем пропорциональной матрице жесткостей. Коэффициент пропорциональности а выбираем из условия равенства логарифмического декремента затухания 0,2 [3].

В нормальной форме Коши система (9) имеет следующий вид:

Ч (0 = Ая (0 + А^ (, я (¿), я (Г)), (10)

где я = {х, х}т — вектор фазовых переменных; А — матрица линейного дифференциального оператора; Ах — вспомогательная матрица; Т () — вектор правых частей системы.

Здесь матрица дифференциального оператора

Л:

Oq Ч

-M-1K -ZM-1D

(11)

где Iу — единичная матрица размерности у х у ; О у — нулевая матрица размерности у х у ; п — число фазовых переменных;

вспомогательная матрица

Ai =

O„ On

On M

(12)

вектор правых частей

f (t)

Jf (t) + R (x, x )J

(13)

где 0 у — нулевой вектор размерности у.

0

n

ч

Общее решение системы (10)

Ч (О = Р (; Яо, ^); Ч ( ) = Ч0, (14)

где ч0 — фазовый вектор в начальный момент времени; t0 — начальный момент времени.

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

решения на одном периоде вектор ч0 полагаем неизвестным и накладываем условие, чтобы все фазовые переменные через период принимали начальные значения. Математически это запишем так:

Ч(Т + ^ = Ч0, (15)

где Т = 2п/ со — период внешнего возмущения. В дальнейшем положим t0 = 0.

С учетом (14) выражение (15) выглядит следующим образом:

Ч0 = р (Т; Ч0,0). (16)

Для решения уравнения (16) используем метод Ньютона [13]: определяем нулевой корень функции

Г (Ч0 ) = р (Т; Ч0,0)-Ч0 = 0. (17)

Параметры рассматриваемой динамической системы, принятые при расчетах, приведены в табл. 1.

Таблица 1

Физические параметры динамической системы

Величина Значение Величина Значение Величина Значение

о, рад/с 314 khx, Н/м 29-108 LG, м 0,030

mW, кг 90,2 кНу, Н/м 27-108 Ар, Па 2 • 106

JEW, кг • м2 1,18 mhx, кг 114,2 S0, мм 0,25

JAW, кг • м2 1,93 mHy , кг 106,3 L1g , м 0,140

mW, кг 90,2 cbl , Н/м3/2 4,4-1011 L 2g, м 0,135

JEW, кг • м2 1,18 Cbl , Н/м3/2 6,7 -1011 JU, Па • с 0,00101

т 2 JAW, кг • м 1,93 Dg , м 0,290

Траектория центра вращения рабочего колеса с учетом влияния ЩУ и без него показана на рис. 9, на котором видно, что ЩУ значительно снижает как статический прогиб вала, так и амплитуду поперечных колебаний рабочего колеса.

Рис. 9. Траектория рабочего колеса

В случае исполнения ЩУ без погрешностей, оно также существенно снижает реакцию в подшипниках, что видно по характеристикам на рис. 10. Нагрузка на исследуемые в данной работе ПК С 2220 К и 22220 ЕК оказывается меньше минимально допустимой, установленной производителем.

Рис. 10. Реакция в правом подшипнике

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

Таблица 2

Варианты расчетов при наличии погрешностей

Номер варианта Параметр

rW , мм Ä, рад F /F0

0 0 0 1

1 0,025 1,2

2 0,130 1,2

3 0,025 16,0е-4 1,2

4 0,250 1,9е-4 2,5

5 0,100 12,0е^ 2,5

Траектории центра вращения рабочего колеса (рис. 11) показаны для каждого из вариантов, приведенных в табл. 2. На рис. 11 видно, что при больших значениях хотя бы одной технологической погрешности рабочего колеса амплитуда колебаний возрастает в десятки раз. Однако даже при малых значениях погрешностей (вариант 1) амплитуда увеличивается более чем в три раза. Следовательно, при исследовании динамики ротора с учетом ЩУ в модель следует включать технологические погрешности.

■200 -150 -100 -50 0 50

V [мкм]

Рис. 11. Траектория центра вращения рабочего колеса с учетом погрешностей

Реакция и момент в ЩУ существенно возрастают (рис. 12) при наличии технологических погрешностей. Реакция в подшипниках (рис. 13) значительно изменяется только при больших значениях погрешностей.

— 6000

CD Q1

0.01 t, [С]

а

б

Рис. 12. Реакция (а) и момент (б) в щелевом уплотнении

Рис. 13. Реакции в правом подшипнике при наличии технологических погрешностей

Спектральные представления характеристик решения приведены на рис. 14. Для случаев 3 и 5, соответствующих большим начальным перекосам рабочего колеса, отмечаются существенные составляющие на удвоенной частоте собственного вращения ротора.

Рис. 14. Спектральные представления характеристик решения

Для насосных агрегатов типа НМ максимальное значение вертикальной составляющей виброскорости vHL опоры, отмеченной стрелкой на рис. 2, составляет 2,8 мм/с по нормам [14]. Среднеквад-ратические значения (СКЗ) vHL для различных вариантов расчетов (см. табл. 2) приведены ниже:

Вариант........................... 0 1 2 3 4 5

СКЗ vHL, мм/с.............. 0,36 0,37 0,25 0,29 0,36 0,41

Для всех вариантов сочетания погрешностей, приведенных в табл. 2, СКЗ vHL существенно ниже значения 2,8 мм/с, что отвечает отличному состоянию насосного агрегата. В то же время амплитуда колебаний рабочего колеса для вариантов 1—3 и 5 значительно возрастает по отношению к случаю, когда погрешности отсутствуют.

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

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

Работа выполнена в рамках проекта «Разработка и производство отечественных насосных агрегатов нового класса для транспорта нефти (импортозамещающие технологии)» при финансовой поддержке Министерства образования и науки РФ (Постановление Правительства РФ № 218).

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

1. Muszynska A. Rotordynamics. - Boca Raton: Taylor & Francis Group, 2005. - 1054 с.

2. Марцинковский В.А. Гидродинамика и прочность центробежных насосов. - М.: Машиностроение, 1970. - 272 с.

3. Бидерман В. Л. Теория механических колебаний: учеб. для вузов. - М.: Высш. шк., 1980. - 408 с.

4. Особенности проектирования упругой компенсирующей муфты: ч. 1. Расчет упругого элемента при соосном расположении валов / Е.О. Билец-кий, А.М. Гуськов, О.А. Ряховский, Е.П. Фирсов // technomag.edu.ru: Наука и образование: электрон. науч.-техн. издание. - 2010. - Вып. 12. -URL: http://technomag.edu.ru/doc/163921.htm

5. Особенности проектирования упругой компенсирующей муфты: ч. 2. Расчет упругого элемента при наличии угловой несоосности валов / Е.О. Билецкий, А.М. Гуськов, О.А. Ряховский, Е.П. Фирсов //

technomag.edu.ru: Наука и образование: электрон. науч.-техн. издание. -2010. - Вып. 12. - URL: http:/http://technomag.edu.ru/doc/163966.html

6. Ковалев М.П., Народецкий М.З. Расчет высокоточных шарикоподшипников. - М.: Машиностроение, 1975. - 280 с.

7. Brennen C. E. Hydrodynamics of pumps. - Norwich: Concepts ETI, Inc.; Oxford: Oxford University Press, 1994. - 293 с.

8. Тараненко П.А. Динамика ротора турбокомпрессора на подшипниках скольжения с плавающими втулками: дис... канд. техн. наук: 05.02.02. -Челябинск, 2011. - 172 с.

9. Itskov M. Tensor algebra and tensor analysis for engineers. - 2d ed. - Dordrecht: Springer, 2009. - 247 с.

10. Ломакин В.О., Петров А.И., Степанюк А.И. Оптимизация геометрических параметров отвода нефтяного магистрального насоса типа НМ // technomag.edu.ru: Наука и образование: электрон. науч.-техн. издание. - 2010. - Вып. 3. - URL: http://technomag.edu.ru/doc/347727.html

11. Создание параметризованных 3D-моделей проточной части центробежных насосов / В.О. Ломакин, П.В. Щербачев, О.И. Тарасов, П.А. Покровский, С.Е. Семенов, А.И. Петров // technomag.edu.ru: Наука и образование: электрон. науч.-техн. издание. - 2010. - Вып. 4. - URL: http://technomag.edu.ru/doc/354657.html

12. Петров А.И., Ломакин В.О. Численное моделирование проточных частей макетов насосов и верификация результатов моделирования путем сравнения экспериментально полученных величин с расчетными // technomag.edu.ru: Наука и образование: электрон. науч.-техн. издание. -2010. - Вып. 5. - URL: http://technomag.edu.ru/doc/356070.html

13. Волков Е.А. Численные методы: учеб. пособие. - 2-е изд., исправл. -М.: Наука. Гл. ред. физ.-мат. лит., 1987. - 248 с.

14. РД 153-39ТН-008-96. Руководство по организации эксплуатации и технологии технического обслуживания и ремонта оборудования и сооружений нефтеперекачивающих станций. - Разраб. 01.01.1997. - Уфа: Ин-т проблем транспорта энергоресурсов, 1997. - 39 с.

Статья поступила в редакцию 28.09.2012