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

замечаем, что (4) удовлетворяется тождественно, а уравнение (3) приводится к виду

д V ,-т2

dz2

+ 2 U 2 V = 0,

где и2 = D2 +ю1 - оператор Даламбера. Отсюда при поперечных колебаниях плиты

V = зт(л/2 П z)F1(х, у); U = ¡ип(>/2 Ц z)

дУ

V = - яп(Т2 Ц z)

дx

1 — — дF

Т х; = л/2 П С08(Л/2 П z)

С1 ду

Из краевых условий (1) следует: — дF

С08(л/2 и Ь)—1 = 0. (27)

дУ

Пусть F:(x,y) - решение уравнения Даламбера

и2 F -4F: = 0.

Тогда, вычисляя (27) при помощи (28), получим

cos 8, = 0.

1 ду

Следовательно, функция сдвиговых напряжений будет отлична от нуля, если cos 81 = 0. п

Значения 81l = —(1 + 2l) (l = 1,2,3,...) представляют собой собственные значения, а функции сдвиговых напряжений F1l(x, у) являются собственными функциями.

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

(28)

Литература

1. Лехницкий С.Г. Анизотропные пластинки. - М.: ГИТЛ, 1947. - 365 с.

2. Лехницкий С.Г. О некоторых вопросах, связанных с теорией изгиба тонких плит // ПММ. -1938. -Т. 2, № 2.

3. Лурье А.И. Пространственные задачи теории упругости. -М.: ГИТЛ, 1955. - 492 с.

Ростовский государственный строительный университет

19 сентября 2006 г.

УДК 621.822+62.251

КОНЕЧНО-ЭЛЕМЕНТНЫЙ АНАЛИЗ ДИНАМИЧЕСКИХ ХАРАКТЕРИСТИК РОТОРНЫХ СИСТЕМ С ПОДШИПНИКАМИ ЖИДКОСТНОГО ТРЕНИЯ

© 2007 г. О.В. Соломин, С.В. Майоров

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

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

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

щихся, параметрических и хаотических колебаний [2 - 4]; 3) способны варьировать в широком диапазоне свои динамические характеристики в зависимости от сочетания рабочих и геометрических параметров [2 - 6].

Поэтому становится очевидной актуальность поиска рациональных параметров роторной системы на стадии проектировочных расчетов. Решение этой задачи на современном этапе не представляется возможным без применения современных САЭ/САЕ-систем. В инженерной практике создания реальных машин, как правило, нет необходимости решения задачи полного моделирования динамики роторной системы. Обычно ограничиваются решением следующих основных задач [6, 7]: 1) определение собственных частот и форм колебаний (модальный анализ); 2) построение амплитудно-частотных характеристик АЧХ в области рабочих частот (гармонический анализ); 3) анализ переходных процессов, позволяющий моделировать импульсные нагрузки на систему (отрыв лопатки, ударное возмущение и т. д.).

Колебания ротора как линейно-упругого изотропного тела, описываются динамическими уравнениями Ламе при соответствующих граничных условиях, которые с учетом диссипативных сил имеют вид [8]

d 2u т E dt2 2 (1 + v)

V 2u +

E

2 (1 + V)(1 - 2v)

grad divu +

„ 2 du u , du

+ 2uV--\— grad div—

^ dt 3S dt

u| r = R-

(1)

где р - плотность; Е - модуль Юнга; V - коэффициент Пуассона; ц - коэффициент внутреннего трения; и - вектор перемещений; / - вектор объемных сил; t - время; Я - заданная функция перемещений на поверхности ротора Г.

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

Граничные условия модели (1) представлены внешними нагрузками и реакциями смазочного слоя подшипников. Внешние силы достаточно легко можно выразить через перемещения, используя закон Гука. Определение реакций смазочного слоя подшипников, строго говоря, требует решения уравнений гидродинамики с учетом значений кинематических параметров цапфы ротора [3, 6, 7, 9]. Однако как показано в ряде работ [2 - 4, 6, 7, 10], учитывая высокие скорости вращения ротора и малые рабочие эксцентриситеты положения центра цапфы ротора, можно использовать линейное представление реакций смазочного слоя подшипников жидкостного трения.

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

X >

Втулка подшипника

Точка 4 Точка 5

Точка 1 ' очка 2 Точка 3

б)

Рис. 1. Динамическая модель опорного узла (а) и модель ротора (б)

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

Rx = Rxо - Kxx AX - Kxy AT - Bxx AX - Bxy AT ; Ry = RT0 - KTX AX - KTT AT - BTX AX - BTT ATT ,

где коэффициенты жесткости и демпфирования подшипника равны:

KXX = -1

KYX = -

BXX = -

dRx dX dRT ~dX

dR

X

dX

KXY = -

KYY = -

BXY = -

dRx dY dRT dY

dRx dY-

Bx =-|dRT

Byy =-drt

d7

тх ~ | I > тт

) о V д ) о

Для решения практических задач динамического анализа удобно рассматривать представление реакций подшипника в форме (см. рис. 1 б):

Ях = Яхо - Сх АХ - Вх АХ ;

Ят = Ят0 - Ст АТ - Вт АТ,

о

о

где соответствующие динамические коэффициенты имеют вид

СХ = KXX + KXY ; BX = BXX + BXY ;

Cy — Kyx + Kyy ; By — Byx + Byy *

Расчет динамических коэффициентов (коэффициентов жесткости и демпфирования) смазочного слоя представляет собой самостоятельную задачу. Определение динамических коэффициентов для конкретного типа подшипника при заданных условиях работы требует, за исключением самых простых случаев, разработки достаточно сложной модели опоры, привлечения численных методов и создания программного обеспечения [3, 6, 7].

Методика численного расчета динамических коэффициентов для различных типов подшипников при турбулентном неизотермическом течении сжимаемого смазочного материала, а также некоторые результаты расчетов подробно изложены в работе [10]. Для решения этой и других задач динамического анализа роторных систем с подшипниками жидкостного трения авторами был разработан программный комплекс АнРоС [11], с помощью которого и были рассчитаны необходимые динамические коэффициенты.

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

Аналитический расчет динамических характеристик такой системы требует решения уравнений с соответствующими граничными условиями типа (1), что является весьма сложной (а в общем случае и невыполнимой) задачей. Применение численных методов делает эту задачу вполне решаемой. На сегодняшний день наиболее популярным методом решения задач динамики реальных конструкций является метод конечных элементов (МКЭ), реализованный в многочисленных программных продуктах (Nastran, Ansys, Cosmos/M и др.). Авторами использовался пакет Ansys, что не снижает общности рассматриваемой методики динамического анализа роторных систем.

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

Прежде всего, при препроцессорной подготовке осуществляется построение трехмерной (3-D Solid) модели ротора. Несмотря на то, что пакет Ansys имеет свой внутренний графический модуль для проведения

трехмерного моделирования, построение сложных реалистичных моделей в нем затруднительно. Поэтому для таких целей предпочтительнее использовать специализированные CAD пакеты (AutoCAD, T-Flex, КОМПАС и др.).

В данном случае в качестве такого CAD пакета был выбран AutoCAD. Моделирование в AutoCAD, как и в любой другой CAD-системе основано на поэтапном моделировании частей ротора (тело ротора, лопатки, шлицы) и дальнейшей их сборке. На рис. 2 а представлены модель ротора и модель лопатки, построенные в AutoCAD (модель лопатки для наглядности увеличена). После создания геометрической модели ротора необходимо ее экспортировать в пакет Ansys, используя, например, формат передачи данных *.SAT.

а)

Рис. 2. Модель ротора: а - 3-0 Solid-модель; б — конечно-элементная модель

После передачи модели в пакет Ansys создается конечно-элементная сетка и задаются граничные условия. Для проведения анализа были взяты следующие элементы (в обозначениях Ansys): 1) для моделирования ротора - <^ОЬЮ92»; 2) для моделирования упруго-демпферных опор - «СОМВШ14».

На рис. 3 а указана предельная форма конечного элемента (тетраэдра) <^ОЬЮ92» с десятью узлами. Каждому узлу соответствует три степени свободы (перемещения по трем координатам). Механические свойства материала этого конечного элемента описываются тремя константами: модулем упругости (Е, Па), коэффициентом Пуассона (V), плотностью (р, кг/м3). Использование данного элемента сочетает приемлемую точность и скорость вычислений. Аппроксимирующая функция для этого элемента имеет вид [12]

и = и1 (Щ —1)Ь 1 + uJ (2Ь2 1)Ь2 ++ик

(2Ьз -1)Ьз +

+и^ (2Ь4 —1)Ь4 + 4(им Ь1Ь 2

+ и^Ь 2 Ь з + ио^1 Ь з +

+u pL~i L 4 + UqL 2 L 4 + uRL 3 L 4

(2)

где и1 - перемещения в узлах; i=(I, ..., К) - обозначения узлов; Ц - нормированные координаты, принимающие значение «0» в узле и «1» на плоскости противоположной узлу; 7=(1, ..., 4) - номер координаты. На рис. 3 а указаны номера плоскостей соответствующие номерам координат Ц. Применяя выбранный элемент <^ОЬЮ92» и выполняя разбиение, получаем конечно-элементную модель ротора на упруго-демпферных опорах (см. рис. 2 б).

Рис. 3. Конечные элементы: а - SOLID92;

б - сомвми

Схема элемента «СОМВЖ14» с двумя узлами изображена на рис. 3 б. Данный тип элемента обладает 12 степенями свободы (по три перемещения и по три угла поворота в каждом узле) и не имеет массы. На рис. 3 б обозначено: k - коэффициент жесткости, Н/м; С„ - коэффициент демпфирования, Н-с/м. Аппроксимирующая функция для этого типа элемента имеет вид

= 2 (uI (1 - ^ )+ uj (1 + 5 ))

(3)

где 5 - локальная координата, принимающая значение «- Ях = Ях 0 - Кхх Ах - Кхт АТ - Вхх АХ - Вхт АТ 1» в узле I и «+1» в узле 3; щ и и3 - перемещения в узлах I и 3 соответственно.

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

Одной из важнейших задач при выполнении конечно-элементного анализа на этапе препроцессорной подготовки является выбор используемых типов конечных элементов и их числа в модели. Определение этих параметров осуществляется на основе согласования двух критериев: требуемой достоверности полученных результатов и времени расчета. Оптимальное

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

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

В качестве параметра верификации был взят прогиб ротора иу в точке приложения вертикальной сосредоточенной силы (см. рис. 1 б, точка 1). Удовлетворительные результаты получаются при использовании в разбиении уже 1000 элементов типа <^01Ю92» (рис. 4, где N - число элементов).

. и, мм

0,03 0,02 0,01 0

х

х

X

I XX

1

10

100 1 103 1 104 1105 1106

N

Рис. 4. Зависимость статического прогиба от числа узлов сетки

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

[м ]^ М+ [В М+ [К Щ = }, (4)

ж2 &

где [М], [В], [К] - матрицы масс, демпфирования и жесткости системы «ротор - подшипники жидкостного трения», определяемые типом и числом выбранных элементов; (2} - вектор перемещений, - вектор нагрузок.

Для каждого типа анализа на основе уравнения (4) получаются разрешающие уравнения. В частности, для модального анализа предполагается, что на систему не действуют внешние силы {Б} = 0. Кроме этого, считаем, что перемещения изменяются по гармоническому закону (2} = (и }е АГ.

Тогда система уравнений (4) может быть записана как следующая однородная система линейных алгебраических уравнений относительно вектора {и}:

n*EQ=3. 007

(5)

(А2 [М] + А[В] + [К])(и} = 0.

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

(2 [м ]+х [в]+[к ])= о. (6)

Заметим, что на практике нет необходимости раскрывать определитель и искать корни полинома. Для отыскания корней уравнения (6) разработан ряд эффективных численных методов [13], которые реализованы в Лтуя. В общем случае, корни алгебраического уравнения (6) А,- - комплексные величины. Значения собственных частот ю,- роторной системы определяют-

ся как

ю г = Im (А г).

Узловые перемещения, соответствующие ,-й собственной частоте можно получить из уравнения (5), подставляя вместо А значение корня А,- уравнения (6) и, положив в векторе {и} любой (например, первый) элемент равным любому вещественному числу (например, единице). По известным значениям собственных векторов {и},- с помощью зависимостей (2) и (3) можно получить форму колебаний (моду), соответствующую 1-й собственной частоте. Заметим, что для обеспечения однозначности определения собственных векторов обычно применяют процедуру нормирования по массе [6]: (и} Т [М](и}, = 1.

При исследовании форм колебаний также часто рассматривается вопрос об устойчивости моды. Считают, что мода устойчива, если выполняется следующее условие (иначе мода неустойчива): р, = Яе(А,)< 0.

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

Построение амплитудно-частотной характеристики (гармонический анализ) осуществляется на основе разделения переменных (координат и времени), при условии воздействия на систему внешних гармонических сил. Предположим, что к системе приложена в некоторой точке сосредоточенная сила Р(Г) = Р^ш(юt + ф), тогда вектор сил можно записать в виде

(^} = (^0 }яп (юГ +ф).

Или, что то же самое, в естественной записи:

а)

б)

в)

{F} = {Fq}6

+ф

г)

Рис. 5. Формы колебаний ротора: а - первая; б - вторая; в - третья; г - четвертая

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

{Q} = {U}вш+ф .

(7)

Подставляя выражение (7) в уравнение (4), получаем следующее уравнение:

(2 [м] + tm[ß]+[K]){U} = {F}.

(8)

Решение уравнения (8) относительно вектора неизвестных {П} выполнятся при дискретном изменении частоты ю в интересующем частотном диапазоне. Затем по полученному массиву векторов {П} достаточно легко построить амплитудно-частотную характеристику для любой точки системы.

На рис. 6 представлены АЧХ точек 2 и 4 (см. рис. 1 б) роторной системы в диапазоне рабочих частот. Эти точки соответствуют вертикальным перемещениям ротора в месте расположения подшипников. По оси абсцисс отложена не круговая частота ю, полученная в Ату&\ а частота колебаний ротора / в герцах в соответствии с выражением ю = 2п(. Как видно из рис. 6,

в области диапазона рабочих частот ротора (до 500 Гц) резонансные зоны отсутствуют, что является позитивным результатом проектирования с точки зрения выполнения критерия вибрационной надежности.

и, мм х 10-3 8,0

7,2 6,4 5,6 4,8 4,0 3,2 2,4 1,6 0,8 0

К, n J , / -—и L.

200

400

600

800 / Гц

а)

U, мм х 10

1,125 1,000 0,875 0,750 0,625 0,500 0,375 0,250 0,125 0

к.

200

400

600

800

/, Гц

б)

Рис. 6. Амплитудно-частотные характеристики: а - для точки 2; б - для точки 4

Анализ переходных процессов проводится на основе интегрирования во временной области уравнения (4) с заданными начальными условиями по методу Ньюмарка [12, 13]. В точке 1 на турбине (см. рис. 1 б) прикладывался ступенчатый импульс силы заданной длительности и величины, и рассчитывались свободные колебания ротора при различных демпфирующих свойствах смазочного слоя. На рис. 7 приведены временные реализации свободных колебаний точки 1 (см. рис. 1 б) при различных значениях коэффициентов демпфирования подшипников, смазываемых жидким водородом. Подробное описание параметров подшипников приведено в работе [10].

U, мм х 10

1,0 0,8 0,6 0,4 0,2 0 -0,2 -0,4 -0,6 -0,8 -1,0

а)

U, мм х 10

2,4 1,6 0,8 0 -0,8 -1,6 -2,4 -3,2 -4,0 -4,8 -5,6

б)

Рис. 7. Отклик на ударную нагрузку для точки 1

На рис. 7 представлены временные реализации колебаний точки 1 при различных давлениях подачи смазочного материала: а) давление подачи номинальное; б) давление подачи - увеличено. В частности, показано [10], что незначительное увеличение давления подачи ведет к малому изменению коэффициентов жесткости, но может увеличить почти в два раза величины коэффициентов демпфирования. Как видно варьирование давления подачи смазочного материала может значительно повлиять на динамические харак-

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

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

Литература

1. Handbook of turbomachinery // - NY, Marcel Dekker, Inc. -

1995. - 472 p.

2. Равикович Ю.А. Конструкции и проектирование подшипников скольжения агрегатов ДЛА. - М.: Изд-во МАИ, 1995. - 58 с.

3. Артеменко Н.П. и др. Гидростатические опоры роторов

быстроходных машин. - Харьков: Основа, 1992. - 198 с.

4. Луканенко В.Г. Колебания высокоскоростных роторов на гидростатических подшипниках и методы снижения виброактивности машин. - Самара: Самарский НЦ РАН, 2001. - 122 с.

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

6. Adams M.L. Rotating machinery vibration: from analysis to troubleshooting. - NY: Marcel Dekker, Inc., 2001. - 354 p.

7. Yamamoto T., Ishida Y. Linear and nonlinear rotordynamics. A modern treatment with applications. - New York, John Willey&Sons, 2001. - 326 p.

8. Вибрации в технике: справочник. В 6 т. - М.: Машиностроение, 1978. Т. 1. - 352 с.

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

10. Соломин О.В. Динамические характеристики гидроста-тодинамических опор в условиях двухфазного состояния смазочного материала // Известия вузов. Машиностроение, 2006, № 1. - С. 14 - 23.

11. Соломин О. В. и др. Анализ роторных систем - АнРоС / Свидетельство об официальной регистрации программы для ЭВМ № 2006610287.

12. Zienkiewich O., Taylor R. The finite element method. Vol. 1. The basis. - Oxford, Butterworth-Heinemann, 2000. - 694 p.

13. Бате К., Вилсон Е. Численные методы анализа и МКЭ. - М.: Стройиздат, 1982. - 448 с.

20 июня 2006 г.

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