CC BY
46
Механика
Вестник Нижегородского университета им. Н.И. Лобачевского, 2012, № 4 (1), с. 171-176
УДК 534.01
МОДЕЛЬ ДИНАМИКИ ГИБКОГО НЕОДНОРОДНОГО РОТОРА НА ЭЛЕКТРОМАГНИТНЫХ ПОДШИПНИКАХ
© 2012 г. В. Ф. Овчинников 1, М.Я. Николаев 1, Е.В. Кирюшина 1, А.А. Кирюшин 1,
В.Н. Литвинов 1, Е.В. Фадеева 1, А.С. Чистов 1, Ф.М. Митенков 2, Н.Г. Кодочигов2
1 НИИ механики Нижегородского госуниверситета им. Н.И. Лобачевского 2 ОАО «ОКБМ «Африкантов», Н. Новгород
Пеступиеа в реСакцию 01.12.2011
Приведена математическая модель, предназначенная для расчетного исследования динамики гибкого неоднородного ротора на электромагнитных подшипниках.
Кеючевые сеева: гибкий неоднородный ротор, математическая модель, система управления, электромагнитный подшипник.
Одним из направлений развития систем подвеса, использующих электромагнитные подшипники (ЭМП), является их применение для гибких роторов сложных конструкций больших габаритов и веса (см., например, [1]). Такие многотонные гибкие неоднородные роторы удовлетворяют сформулированным в работе [2] критериям, определяющим сложную уникальную систему. В подобных роторах соединены разнородные машины: турбины, электрогенераторы, компрессоры и т.д. Каждый ЭМП имеет свою систему управления, взаимодействующую с системами других подшипников через ротор. Многочисленные силы, действующие на ротор (вес, силы от дисбаланса, гироскопические, термомеханические, газодинамические, электромагнитные силы), возбуждают одновременно разные формы колебаний гибкого ротора, комбинация которых может привести к существенно различающемуся характеру движений разных частей ротора. Система в целом сугубо нелинейная. Для исследований, необходимых как при проектировании, так и при эксплуатации подобных систем, построена модель динамики гибкого неоднородного ротора на ЭМП, которая включает: механическую модель, модели сил разной физической природы, модель системы управления.
Механическая модель базируется на уравнениях и результатах исследований динамики гибких роторов [3]. В большинстве задач динамики ротор рассматривается как упругий неоднородный стержень. При его движении можно выделить четыре вида колебаний: как твердого тела, упругие крутильные, упругие продольные и упругие изгибные колебания в двух направлениях. В области малых колебаний для прямых
стержнеи упругие крутильные, продольные и изгибные колебания независимы. Особенности динамики ротора на ЭМП обусловлены изгиб-ными колебаниями, поэтому эти движения рассматриваются в работе.
Для формирования механической модели ротора используется прямоугольная система координат 0ху,. Ось 0х направлена вертикально и совпадает с осью вращения недеформирован-ного ротора. Ротор разбивается на серию соединенных между собой однородных участков кругового сечения. Для описания изгибных колебаний участка используется модель балки Тимошенко [4], в соответствии с которой определяющие динамику ротора функции связаны уравнениями
дВ, д’-иг дм, ,д20,
"аГ=т1Г ~ ч’ • ~зГ =~в,+р'~дг
1
— , ---- = © +—-Т В .
дхЕІ дх 2 Gfr y
а© 1
-=— м.
(1)
Функции Ц,(х,0, 02(х,О характеризуют линейные и угловые перемещения ротора при движении в направлении оси 0у, функции Мг(х,0, Qy(x,t) - возникающие в сечении ротора внутренние моменты и перерезывающие усилия. цу, - действующие на ротор распреде-
ленная сила и распределенный момент; Е, V, С=Е/2(1-у) - модуль упругости, коэффициент Пуассона и модуль сдвига; р - плотность материала ротора; f, I - площадь и осевой момент инерции нормального сечения участка; т - погонная масса участка; г* - коэффициент, учитывающий неравномерность распределения касательных напряжений по нормальному сечению участка [5]. Для объединения в механиче-
скую модель используется два типа соединения отдельных участков: жесткое и опора. Функции иу(х,0, ©г(х,0 и М2(х,?) в сечениях соединения участков непрерывны. Функция Qy(x,t) в сечении жесткого соединения участков непрерывна, а в сечении расположения опоры терпит разрыв, величина которого равна силе реакции опоры.
Уравнения (1) определяют изгибные колебания ротора в направлении оси 0у, аналогичные уравнения можно записать и для колебаний ротора в направлении оси 02, которые характеризуются функциями Ы2(х,(), ©у(х^), Му(х,(), Qz(x,t).
Для перехода к дискретной модели ротора определяющие динамику ротора функции разлагаются по базовой системе функций - формам собственных колебаний Ц(х), ©к(х), Мкх), Qk(x). Для перемещений это разложение имеет вид:
иу =^ак(tЫ(х), и2 = ^Ьк(t)ик(х). (2) к=1 к=1
К - количество форм собственных колебаний в аппроксимации решения. Аналогичные выражения используются для функций углов поворота, моментов и перерезывающих усилий.
Выбор базовой системы функций определяется решаемой задачей. Для расчета динамики упругих конструкций при наличии нелинейных внешних связей можно использовать формы собственных колебаний безопорного ротора в комбинации с формами движения ротора как твердого тела. В работе [6] для моделирования гибкого неоднородного ротора на ЭМП использовались формы собственных колебаний ротора с шарнирным опиранием в местах расположения ЭМП в комбинации с формами статических деформаций. В данном случае в системе управления ЭМП предполагается формировать требуемые законы (в том числе линейные) изменения сил в зависимости от перемещений ротора [7, 8], поэтому в качестве базовой системы функций целесообразно использовать ортогональные формы собственных колебаний ротора при граничных условиях, соответствующих опиранию ротора на упругие опоры жесткости сп в сечениях х = хп расположения ЭМП и свободных концевых сечениях (п = 1, N, N - количество радиальных ЭМП). Условия ортогональности и нормировки базовой системы функций имеют вид
}(тики1 + р!®к®J1 = Р кк7=7.|; к,] = \К, (3)
где т0 - масса ротора.
Интегрирование выполняется по всей длине ротора I.
Коэффициенты разложения в выражениях (2) играют в механической модели ротора роль
обобщенных координат. Для получения разрешающих уравнений используются уравнения Лагранжа второго рода [9].
Кинетическая энергия колебаний и потенциальная энергия деформаций при движении ротора в направлении оси 0у определяются выражениями
После подстановки полученных выражений в уравнения Лагранжа с учетом свойств форм собственных колебаний математическая модель динамики ротора в направлении 0у примет вид
1 а 0 /.л
т0 + т°а а = Ка . (5)
В этой записи а, Ь - векторы обобщенных координат размерности К, элементами которых являются коэффициенты разложения (2),
О = diag (ю2,ю2,...,юК)- диагональная матрица порядка К, элементами которой являются квадраты собственных частот ротора.
Компоненты векторов обобщенных сил определяются стандартным образом [9], и для них с учетом уравнений (1) и разложения (2) можно записать
I
Ка = {#} К =\^уик(х) + Цг©к(х)К к = 1К.
0
Формально эта запись справедлива и для сосредоточенных воздействий, которые можно описывать с использованием дельта-функций. Ниже представлены выражения для наиболее значимых сил, учтенных в математической модели.
В заданном положении ротор удерживают силы Fn радиальных ЭМП, которые зависят от горизонтальных перемещений ротора йп = = иу(хп) в сечениях расположения ЭМП и от токов в катушках. В силу (2) для вектора перемещений й = {йп} справедливо представление 1 н = {}гпк ^ ^ = ик (хп).
Знание вектора перемещений и токов в катушках ЭМП дает возможность определить вектор сил F = {Fn}. Этот вектор можно представить в виде суммы
F =-ей + ДF, с = diag (с1,..., ^).
Первая составляющая этой силы характеризует реакцию упругих опор, и она уже учтена в выражении (4) для потенциальной энергии, и следовательно, и в уравнениях (5). В связи с этим в уравнениях (5) следует учесть только
обобщенную силу, порождаемую вторым слагаемым выражения
Ra = НТ ЛР = Са + НТР, С = НтсН.
Верхний индекс «Т» соответствует операции транспонирования матрицы.
Гироскопические силы в уравнениях движения (1) учтены распределенными моментами
д©, д©
М, = -2юр1—, м , = 2юр1 ~~, д, = д, = °
ді
ді
dф
ю =-------угловая скорость вращения ротора, ф
dt
- угол поворота ротора. С учетом разложения (2) для векторов обобщенных сил получены выражения:
к, ] = 1, К.
Наличие дисбаланса соответствует действию на ротор распределенной силы
д, = ю2т[е1( x)cos ф- e2(x)sin ф], д, =aгm\el(x)sin ф + e2(x)cos ф], ц, = ц = 0.
Функции еі(х), е2(х) характеризуют положение центра масс нормального сечения ротора по отношению к оси вращения в системе координат, связанной с ротором. Для соответствующих векторов обобщенных сил получены выражения:
= Ю^^Ф - 52sinф], 5 = К,},
I
Які = ^теі (х)ик (х)сХ, к = 1, К, і = 1,2.
0
По аналогичной схеме получены выражения для обобщенных сил другой физической природы. Для вертикально расположенного ротора влияние силы тяжести на поперечные движения ротора в уравнениях движения учитывается распределенными силами:
дх
д00© ,
дх
М2 = М у = 0.
В этом выражении Q0(x) - действующая на ротор осевая статическая сила.
Электромагнитные силы в генераторе, возбудителе и двигателе, а также газодинамические силы в турбине, компрессорах и в лабиринтных уплотнениях их валов моделируются силами тяжения
д = г и , д = г и , и = и = 0
“у е у’ е г? г^у
и циркуляционными силами
д = г и , д =—ги , и = и = 0,
*2 у V г? V у? г^у ’
ге, Гу, - коэффициенты, зависящие от конструкции и характеристик указанных элементов.
Одним из основных факторов, существенно влияющих на динамику ротора, является диссипация энергии, обусловленная внешним и внутренним трением [3]. Силы внешнего вязкого трения, связанные с взаимодействием ротора с окружающим газом, зависят от конструкции ротора и его скорости. Аналогичные по своему демпфирующему влиянию на колебания ротора электромагнитные силы в ЭМП по величине значительно превосходят силы внешнего трения, поэтому учет сил внешнего трения не приводит к заметному влиянию на точность результатов расчетного анализа динамики ротора.
Внутреннее трение может оказать существенное дестабилизирующее влияние на динамику ротора. Многочисленные исследования показали [10], что в широком диапазоне частот (0-10000 Гц) уровень рассеяния энергии в конструкционных материалах не зависит от скорости деформаций, а определяется видом напряженного состояния, уровнем деформаций, температурой материала. Аналогичными свойствами обладает и конструкционное демпфирование в узлах соединения. В настоящее время имеется большое количество теорий внутреннего трения, которые в отдельных случаях дают возможность учесть основные особенности процесса внутреннего рассеяния энергии. При этом используется связь между тензорами деформаций е и напряжений ст в процессе колебаний (в том числе и нестационарных) ст=ст(е). В общем случае эта зависимость нелинейная и неоднозначная (гистерезисного типа). Универсальные модели такой зависимости отсутствуют. Большинство теорий и моделей внутреннего рассеяния энергии, которые согласуются с экспериментами, очень сложны в практическом применении или ориентированы на определенный круг частных задач (установившиеся колебания, одночастотные колебания, линейные задачи) [11,12]. В связи с этим при моделировании динамики ротора для учета внутреннего трения избран путь компромиссного решения, который основан на двух положениях: а) заметное влияние на процесс колебаний рассеяние энергии оказывает вблизи резонансов при движении ротора на частотах, близких к собственным; б) на характеристики колебательных процессов принципиальное влияние оказывает уровень рассеяния энергии за период колебаний (площадь петли гистерезиса), при этом детальная зависимость напряжений от деформаций (форма петли гистерезиса) имеет второстепенное значение.
В разработанной модели динамики гибкого неоднородного ротора внутреннее трение для всех обобщенных координат учитывается неза-
висимо. Для каждой обобщенной координаты используется феноменологическая модель фрикционного демпфирования (модель Корчин-ского, модель Леонова и Безпалько) [13]. Наличие в материале ротора внутренних потерь описанного вида приводит к появлению в каждом уравнении движения для обобщенной координаты qi, характеризующей отклонения ротора в системе координат, вращающейся с ротором, дополнительной обобщенной силы
R> = -r x| qt (“sign , і = 1, к, (6)
dt
где ri - парциальная жесткость, соответствующая координате q{; х, а - коэффициенты, характеризующие уровень рассеяния энергии и зависимость от амплитуды колебаний.
При циклических деформациях с амплитудой Aqi потери энергии при наличии силы (6) определяются как
4 r X
A W =
-A qa
8X a +1
Aq,a
а +1
Общий уровень колебательной энергии характеризуется величиной !¥= г Ддг2/2, а относительное рассеяние энергии - отношением
(7)
Соотношение (7) с учетом связи логарифмического декремента колебаний 8 с относительным рассеянием у (у = 28) дает возможность по экспериментальным данным для логарифмических декрементов определить параметры х и а. Отметим, что при а =1 соотношение (6) дает модель амплитудно-независимого демпфирования, как и в модели линейного вязкого трения. В этом случае
R. =-г.хф'(діX ф"(ді) = \qt|sign <dq-, x = s/2. (8)
dt
Соотношению (8) соответствует петля гистерезиса, состоящая из треугольников в осях Ri и q..
Векторы обобщенных координат a, Ь (2) характеризуют движение ротора в неподвижной глобальной системе координат Gxyz. При этом вектор a определяет движение ротора в направлении оси Gy, а вектор Ь - в направлении оси Gz. Для описания деформаций ротора следует рассмотреть движение ротора в системе координат GXYZ, жестко связанной с ротором. Ось GX подвижной системы координат совпадает с осью Gx глобальной системы координат, оси GY и GZ образуют угол ф соответственно с осями Gy и Gz глобальной системы координат. Принимая во внимание связь между подвижной и неподвижной системами координат, для векторов обобщенных координат A, B, характеризующих
движение ротора в подвижной системе координат в направлении осей 07 и 02 можно записать: А = а То8ф + Ъ 8тф; В =-а 8Шф + йс То8ф. (9) Процедура учета внутреннего рассеяния энергии сводится к введению в соответствующие уравнения движения диссипативных сил т т
ГА = - -^ДШ>(А), Рв = - -^Ф(В), (10),
где Д - диагональная матрица логарифмических декрементов затухания для форм собственных колебаний. Запись Ф(д) характеризует векторную функцию, компоненты которой определяются соответствующими компонентами вектора д:
ф={ф,}, д=д}, йд={й^}, фі = ф*(д,).
Обобщенные силы (10), записанные в подвижной системе координат, дают добавки в правые части исходных уравнений (5) (обратное преобразование для соотношений (9))
Яа = р ТО8ф - РВ 8Шф, ЯЪ = Р 8ІПф + р ТО8ф. (11) Таким образом, учет внутреннего рассеяния энергии в материале ротора и учет конструкционного демпфирования при деформациях ротора сводится к учету в математической модели динамики ротора (5) обобщенных сил (11).
Сила Рм, действующая на ротор со стороны магнита, определяется как производная по смещению ротора от магнитной энергии, сосредоточенной в зазоре между обмоткой электромагнита и ротором [14]. Сила зависит от смещения ротора Ур и от тока 1м, протекающего в обмотке магнита. В свою очередь, ток зависит от управляющего напряжения им, подаваемого на обмотку магнита ЭМП, а также от ее индуктивности L и активного сопротивления Я:
Г г V
F, = G,
LgSg
2
V So - Yp J
L=
LoS
S - yd
d^M
_____м
dt
+ RI = U
L0 - индуктивность обмотки магнита при центральном расположении ротора; Тм - полный магнитный поток (потокосцепление) в обмотке магнита; Gм - коэффициент, зависящий от конструкции ЭМП, используемых материалов, 5 -эффективный (магнитный) зазор, характеризующий зависимость текущей индуктивности обмотки магнита от смещений ротора; 50 - эффективный зазор, характеризующий зависимость силы, действующей на ротор со стороны обмотки магнита, от смещений ротора.
Задача системы управления - формирование управляющего напряжения, подаваемого на обмотку магнита.
м
А (мм)
0.25 -
А( мм)
0.25 •
° “----------------
/. С—Л
Е(Гц) ~-1---1---1---1-
Г(Гц)
0 40 80 120 160 О 40 SO 1.20 LS0
Рис. 1. Амплитуда перемещений ротора в месте нахо- Рис. 2. Амплитуда перемещений ротора в месте нахож-ждения верхнего ЭМП в режиме с постоянными коэф- дения нижнего ЭМП в режиме с постоянными коэффи-фициентами СУ и в режиме со сменой коэффициентов циентами СУ и в режиме со сменой коэффициентов
СУ при подходе к критическим частотам
СУ при подходе к критическим частотам
Исходной информацией для системы управления являются перемещения ротора у, в месте расположения датчика. Перемещения ротора преобразуются в сигнал перемещений Уе с учетом постоянной времени датчиков Т\.
Т + У = У,.
1 7, е а
dt
Сигнал перемещений задается последовательностью Ук через дискретный период опроса датчиков перемещений Т2:
Ук = Ус(д, ^ ^ + кТ2, к = 0,1,2, ..., (12)
где ^ - сдвиг во времени, учитывающий неод-новременность опроса датчиков разных ЭМП.
Последовательность (12) является исходной информацией для формирования тока рассогласования I* и тока управления 1упр. При использовании пропорционально-интегрально-дифференциального регулятора
1*(1)=Лпр Ук-1 +Вдиф (Ук-1 - Ук-2)/ Т2 +
+ Синт Т2 (Ук-1 + Ук-2 +■■■), tk <^< tk+1,
*упр
- Io +1* или IyV = Io -1*,
Лпр, Вдиф, Синт - пропорциональный, дифференциальный и интегральный параметры регулятора, 10 - номинальный рабочий ток в обмотке магнита.
Ток управления ограничивается по величине:
ССЛи 1упр > 1тах , то 1упр 1тах ;
ССЛи 1упр < ° то 1упр °
1тах - максимальный ток в обмотке электромагнита ЭМП.
Ток управления является исходной информацией для формирования управляющего напряжения им, подаваемого на обмотки магнита. Управляющее напряжение может принимать три значения: 0, +и„ох, -итах.
Алгоритм формирования управляющего напряжения следующий.
Вкеючение, те есть изменение им ет иуея
Се +итах иеи итах:
Есеи Яес(1м - 1упр) < ДU, те им +итах-
Есеи Яес(1м - 1упр) > +ДU, те им итах-
Выкеючение, те есть изменение ет +итах иеи -итах Се нуея:
При им = +итах, есеи Яес(1м - Іупр) > +Ди, те им =0.
При им = -итах, есеи Яес(Ім - Іупр) < -Ди, те
им =0.
Яес - коэффициент обратной связи, имеющий размерность сопротивления; Ди - параметр, характеризующий зону нечувствительности системы управления.
Таким образом, в модели системы управления ЭМП учтены следующие основные факторы, степень влияния которых на динамику ротора на ЭМП подлежит расчетному и экспериментальному исследованию:
- нелинейная зависимость электромагнитных сил от тока в обмотках ЭМП и смещений ротора;
- отличие координат расположения датчиков и соответствующих ЭМП;
- инерция электромагнитов;
- частота съема информации с датчиков положения ротора в пространстве;
- время формирования управляющего напряжения (от момента опроса датчика до поступления управляющего напряжения на обмотки ЭМП);
- зона нечувствительности и гистерезиса системы управления;
- ограничения на максимум тока в обмотках ЭМП.
Пример расчета динамики гибкого вертикального ротора массой т = 10.8 кг на двух радиальных и одном осевом ЭМП приведен на рис. 1, 2. Зазор между ротором и резервным подшипником равен 0.2 мм. На графиках даны
зависимости амплитуды радиальных перемещений ротора в месте расположения верхнего и нижнего ЭМП (А) от частоты вращения (Р) при разгоне при постоянных коэффициентах системы управления Лпр и Вдиф в течение всего разгона (черные кривые) и со сменой коэффициентов Лпр и Вдиф при подходе к критическим частотам вращения 48 Гц, 124 Гц и при частоте 110 Гц (красные кривые).
Работа выполнена при финансовой поддержке Министерства образования и науки РФ и Российского Фонда Фундаментальных исследований (грант 10-0800882 а).
Список литературы
1. Митенков Ф.М., Кодочигов Н.Г. и др. Высокотемпературный газоохлаждаемый энергоисточник для промышленного производства водорода // Атомная энергия. 2004. Т. 97. Вып. 6. С. 432-446.
2. Митенков Ф.М., Знышев В.В. и др. Проблемы и принципы математического моделирования динамики сложных уникальных систем // Математическое моделирование. М.: Наука, 2007. Том 19. № 5. Стр. 39-44.
3. Вибрации в технике. Справочник в 6 томах. Том 3. Колебания машин, конструкций и их элементов / Под ред. Ф.М. Диментберга и К.С. Колесникова. М.: Машиностроение, 1980. 544 с.
4. Тимошенко С.П. Колебания в инженерном деле. М.: Физматгиз, 1950. 436 с.
5. Филиппов А.П. Колебания деформируемых систем. М.: Машиностроение, 1970. 736 с.
6. Знышев В.В., Кирюшина Е.В., Николаев М.Я. и др. Моделирование динамики вертикального неоднородного гибкого ротора на электромагнитном подвесе // Вестник Нижегородского университета им. Н.И. Лобачевского. Серия Механика. 2006. Вып. 1(7). Н. Новгород: Изд-во Нижегородского госуниверси-тета. С. 14-19.
7. Знышев В.В., Кирюшина Е.В., Литвинчук С.Ю. и др. Алгоритм формирования заданной силы электромагнитных подшипников в системе управления электромагнитного нодвеса ротора //Вестник Нижегородского университета им. Н.И. Лобачевского. 2010. Вып. 5(1). С. 138-141.
8. Баландин Д.В., Коган М.М. Управление движением вертикального жесткого ротора, вращающегося в электромагнитных подшипниках // Известия РАН. Теория и системы управления. 2011. № 5. С. 3-17.
9. Гантмахер Ф.Р. Лекции но аналитической механике. М.: Наука, 1966. 300 с.
10. Хаджиян, Марси, Сауд. Обзор методов нахождения оценок эквивалентного демпфирования по данным экспериментов // Труды Американского общества инженеров-механиков. Теоретические основы инженерных расчетов. М.: Мир, 1988. №1. С. 163-175.
11. Писаренко Г.С. Колебания механических систем с учетом несовершенной упругости материала. Киев: Наукова думка, 1970. 380 с.
12. Сорокин Е.С. К теории внутреннего трения нри колебаниях упругих систем. М.: Госстройиздат, 1960. 132 с.
13. Пановко Я.Г. Внутреннее трение нри колебаниях упругих систем. М.: Физматгиз, 1960. 196 с.
14. Schweitzer G., Bleuler H., Traxler A. Active magnetic bearings. Basics, Properities and Applications. Zurich, Vdf Hochschulverlag AG an der ETH, 1994. 244 p.
A MODEL OF THE DYNAMICS OF A FLEXIBLE NONUNIFORM ELECTROMAGNETICALLY SUSPENDED ROTOR
V.F. Ovchinnikov, M.Ya. Nikolaev, E. V. Kiryushina, A.A. Kiryushin,
V.N. Litvinov, E. V. Fadeeva, A.S. Chistov
A mathematical model for numerical simulation of the dynamics of a flexible nonuniform rotor on electromagnetic bearings is presented.
Keywords: flexible nonuniform rotor, mathematical model, control system, electromagnetic bearing.
