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

Наука и Образование

МГТУ им. Н.Э. Баумана

Наука и Образование. МГТУ им. Н.Э. Баумана. Электрон. журн. 2016. № 08. С. 202-233.

]Э5М 1994-040В

Б01: 10.7463/0816.0844628

Представлена в редакцию: Исправлена:

© МГТУ им. Н.Э. Баумана

УДК 62-523; 612.13

Влияние внешнего демпфирования на динамику жёсткого ротора на магнитных подшипниках

04.07.2016 18.07.2016

Скорюков С. В. , Гуськов А. М.

асш£коу_ат @таД.ги 1МГТУ им. Н.Э. Баумана, Москва, Россия

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

Ключевые слова: механическая поддержка кровообращения (МПК), искусственный желудочек сердца (ИЖС), активные магнитные подшипники (АМП), динамика жёсткого ротора, устойчивость движения, бинормальные координаты

1.Введение

Задача о движении тела в сопротивляющейся среде (например, о падении тела в воздухе) интересует исследователей вот уже несколько столетий: ещё в Средние века появилась необходимость изучения зависимости дальности стрельбы от величины угла возвышения ствола пушки [1]. Как показывают эксперименты, прямолинейное стационарное свободное торможение тела в среде неустойчиво. Но также следует учитывать дополнительный параметр, характеризующий вращательную производную момента по угловой скорости тела. Этот параметр вносит в систему диссипацию, увеличивая устойчивость движения обтекаемого тела [1].

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

крови из полости желудочка в восходящий отдел аорты. В клинической практике такой насос называется искусственным желудочком сердца (ИЖС) [2, 2].

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

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

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

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

Существует два основных способа управления АМП: раздельное и многосвязное (оптимальное). В первом способе каждый из каналов управления регулируется автономно по сигналу своего датчика положения, перекрестные связи между каналами отсутствуют. Такой тип управления используется в подавляющем большинстве случаев ввиду простоты регулятора и удобства его наладки. При оптимальном управлении существуют связи между каналами.

АМП находят широкое применение в станках, вакуумных системах и системах очистки помещений, турбомашинах, медицинском оборудовании [5] ввиду их неоспоримых преимуществ [6] по сравнению с подшипниками, традиционно используемыми в технике:

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

• наличие электронной системы управления (контролируемость положения оси, регулируемость жёсткости и демпфирования подвеса).

• АМП обладает также и рядом недостатков:

• необходимость во внешнем источнике энергии;

• относительно низкая несущая способность;

• необходимость в страховочных подшипниках на случай отказа АМП [7];

• сложность электронного блока управления и некоторыми другими.

Из этого вытекает, что АМП находят применение там, где возможности традиционных подшипников качения исчерпаны. Подробнее о преимуществах и недостатках АМП можно узнать в работе [6].

Активные магнитные подшипники - существенно нелинейная система. Нелинейность в АМП обусловлена следующими факторами [5]:

1) нелинейная зависимость магнитной силы от смещения и от силы тока в катушке;

2) насыщение магнитного потока при высоких значениях тока в катушке или закрытии зазора между валом и корпусом подшипника;

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

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

2. Цель исследования

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

3. Исходные данные

Исходные данные были приняты с ориентацией на работу [2].

Исходные данные, характеризующие материал ротора представлены в таблице 1:

Таблица 1. Характеристики материала ротора

Параметр Символ Размерность Значение

Материал ротора - - Титан с магнитными вставками

Плотность Р кг ■ м3 4540

4. Расчётная схема, ограничения

Рассмотрим симметричный однородный жёсткий ротор массой т как стержень, установленный в упруго-демпфирующих опорах (АМП1 и АМП2). Ротор вращается с постоянной угловой скоростью О.

Представление АМП упруго-демпфирующими элементами уместно, если пара датчик-усилитель является совмещённой [5]. В нашем случае рассматривается ротор с несовмещённой парой. Поэтому необходимо учесть несовпадение положения датчика и подшипника.

Модель ротора представлена на рисунке 1:

beanncj plane

sensor plane "¿l"

Ьеачпа plane "a"

sensor plane "c "

Рисунок 1. Неподвижная и связанная системы координат для описания модели жёсткого ротора

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

Для решения поставленной задачи оговорим ряд ограничений [5]:

• ротор симметричный и жёсткий;

• в исходном положении ротора, когда он находится в покое, центр масс совпадает с началом неподвижной системы координат I - х1у121 (см. рисунок 2);

• отклонения от начальных размеров малы по сравнению с размерами ротора - это позволяет линеаризовать уравнения движения и разделить радиальное и осевое перемещение;

• положение ротора, учитывающее поступательное и угловое смещения, характеризуются положением связанной с ротором системы главных осей 5 - х0у020 по отношению к I - х1у121;

• угловая скорость О относительно продольной оси г0 предполагается постоянной;

• движение ротора описывается перемещениями х8, у8 центра масс S по отношению к неподвижной системе отсчета I - ХУ^; отклонения и угловое вращение ротора вокруг оси описываются тремя карданными углами а, Р, у, где П = у.

Линеаризация приводит к двум характерным углам вращения а, Р вокруг осей х1 и у соответственно.

Поскольку ротор жёсткий, система имеет 4 степени свободы: перемещения х, у в каждой из двух опор (см. рисунок ниже). Вектор перемещений примет вид [5,8]:

Ян = {хл хв Ул Ув Г (!)

Рисунок 2. Ввод систем координат для модели жёсткого ротора

Но математически удобнее решать систему в координатах, привязанных к центру тяжести ротора. Таким образом, положение ротора описывается вектором обобщённых координат:

Я = {р х5 -а у5 }Г (2)

Информация о положении центра масс ротора приведена в таблице 2: Таблица 2. Информация о линейных размерах относительно центра масс

Параметр Символ Размерность Значение

Длина 1 м 22 ■1(Г3

Диаметр ротора ¿1 м 15,6-10"3

Расстояния от центра масс до подшипников а Ь м -9 ■ 1(Г3 9■10~3

Расстояния от центра масс до сенсоров с й1 м -10-Ю"3 10-Ю"3

5. Уравнения движения жёсткого ротора, закреплённого в АМП

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

Уравнения, описывающие поведение системы, будем выводить в системе координат, связанных с центром масс (2).

Уравнение движения системы:

Щ + С^ = ВГ (3)

Уравнение движения системы с учётом присутствия сил дисбаланса:

Мч+(541 = 81+^ (4)

где М - симметричная положительно определённая матрица инерции системы:

1у 0 0 0

0 т 0 0 0 0 1х 0 000т

M

(5)

G - кососимметричная матрица гироскопических моментов:

0 0 10"

G = J Q

(6)

B

(7)

0 0 0 0 -10 0 0 0 0 0 0

В — матрица связи вектора сил в привязке к координатам АМП q в и вектора сил в привязке к координатам центра тяжести q :

" a b 0 0 110 0 0 0 a b 0 0 11

fu - вектор внешнего возмущения. Примем синусоидальное возмущение с амплитудой A и частотой QExt:

fu ={0 Af ■ sin (üExtt) 0 0}r (8)

f - вектор нагрузок в системе координат, связанных с АМП:

f = {fxA fxB fyA fyB} (9)

Для дальнейшего исследования необходимо описать природу силы f . Поскольку была поставлена задача учесть несовпадение положений датчиков и подшипника, кроме векторов положения q, qB, введём вектор, характеризующий положение ротора относительно датчиков:

У = {XseA XseB ySeA ySeB } (Ю)

Связь вектора y и q определяется следующим выражением:

y = Cq (11)

где С - матрица связи вектора положения центра тяжести ротора в привязке к системе ко ординат я и вектора положения центра тяжести ротора относительно датчиков у :

" с 10 0

а 100

0 0 с 1

00 а 1

с =

(12)

6. Вязко-упругая интерпретация процессов в АМП и их учёт в уравнении

движения

Как было сказано выше, АМП является существенно нелинейной системой. Рисунок 3 иллюстрирует нелинейную зависимость магнитной силы от величины зазора 5 и токов в катушке электромагнита /' [5].

Рисунок 3. Нелинейные зависимости управления силой

Как показывает практика, проблема учёта нелинейности может быть успешно решена путём линеаризации зависимости магнитной силы от величины смещения ротора от положения равновесия и токов в катушке (рисунок 4). Для случая магнитного подвеса

двухстороннего действия, ток в катушке 1 = {¡ы ¡уВ} [5], тогда выражение для

магнитной силы можно записать в виде:

Г 1 ) = -К>Я в + К11

(13)

K.

KA 0 0 0 " ~ кiA 0 0 0

0 ksB 0 0 , K. = 0 kiB 0 0

0 0 Ka 0 ? 1 0 0 kA 0

0 0 0 ksB _ 0 0 0 kiB

(14)

Матрица Ks (force-displacement factor) связывает магнитную силу и перемещение тела. Как и жесткость пружины, работающей на растяжение-сжатие, ks измеряется в

Н ■ м _ 1 Однако, если сила упругости пружины всегда направлена в сторону, противоположную перемещению, и стремится вернуть тело в положение равновесия, то рассматриваемая составляющая магнитной силы направлена в сторону перемещения, и, следовательно, стремится дестабилизировать систему. Поэтому Ks называют матрицей отрицательных позиционных жесткостей магнитного подвеса ( ks < 0 ). Матрица K. (force-current factor) — матрица коэффициентов к (токовых жесткостей подвеса), измеряемых в Н ■ А _ 1.

Она связывает магнитную силу и ток [5,6].

Упругие и демпфирующие свойства активного магнитного подвеса определяются его конструктивными особенностями и силой подаваемого тока. Согласно работе [2], общем случае отрицательная жёсткость и токовая жёсткость соответственно:

КА,В = -Ls '%4A,B ; kAB = ; S = 1, 2, 3, 4 (15)

где Ls - индуктивность активных обмоток, Гн; ns - коэффициент числа полюсов 5 - зазор, м; i - ток смещения, . Положим далее, что соответствующие коэффициенты жёсткости в опорах равны:

ksA - ksB - kS ; kiA - kiB - ki

(16)

После объединения уравнений (4), (13) и совершая переход в координаты q :

qB=BT q (17)

получим:

Обозначим матрицу коэффициентов отрицательных позиционных жёсткостей в координатах центра тяжести (COG) [5]:

Kss - BKSBT (19)

Поскольку матрица Ks диагональная, матрица KSs является симметричной. Перенося

слагаемое, включающее Ks, в левую часть, получим уравнение, содержащее справа

только управляющее слагаемое, а слева - слагаемые, содержащие симметричные матрицы:

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

7. Децентрализованное управление магнитным подвесом

Обратимся к модели одностороннего подвеса. Поскольку АМП можно представить в виде вязко-упругой модели [5], имеет место выражение:

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

Из этого выражения получим закон изменения тока в зависимости от положения ротора:

Расширяя выражение на случай двустороннего подвеса, имеем:

Это выражение описывает децентрализованное (местное) ПД-управление. Это означает, что подшипники и сенсоры не связаны, то есть их оси разнесены на некоторое расстояние. Такой способ используется в большинстве промышленных систем с магнитными подшипниками, сенсоры обычно не встраиваются в подшипники [2,5]. Схема ПД-управления представлена на рисунке ниже:

Рисунок 5. Схема децентрализованного ПД-управления

Учитывая (11):

1=-РтСч-ОтС,1 (24)

Здесь Рт, Бт матрицы регулировочных коэффициентов или коэффициентов усиления обратных связей по перемещению и по скорости соответственно, которыми описывается PD-управление [5]:

Я 0 0 0

0 Рв 0 0 0 0 ^ 0 0 0 0 Рв

ВА 0 0

0 0

0

0 Бл

0 0 0

0 0 0

(25)

Коэффициенты матрицы Рт измеряются в А ■ м 1 , а коэффициенты матрицы Бт - в

. Приведем ориентировочные формулы для определения коэффициентов управления [2]:

Р = Р =-

1 л 1 в

2к

Бл = Бв =

V

т •К 2

(26)

к/ л в к/

Преобразуем систему уравнений движения ротора для случая PD-управления. Подставим выражение (24) в (20)

Мч + Сц + К^ч = ВК 1(РШСЧ +- О.С,) + Гп (27)

Введём обозначения:

к = вк.Рс

С , ш (28)

Бс = ВХ,БтС

К, соответственно, обеспеченные PD-управлением матрица жёсткости и матрица демпфирования.

Введём обозначения:

К = +

88 с (29)

^ = С + Бс

К - матрица при векторе обобщённых координат q . - матрица при векторе производных вектора обобщённых координат ц, не учитывающая внешнее вязкое демпфирование.

Перенося в левую часть равенства, окончательно имеем уравнение движения ротора в координатах q:

8. Ограничения децентрализованного управления

Как было показано ранее, использование децентрализованного (местного) управления при моделировании динамики жесткого ротора в активных магнитных подшипниках, физически оправдано и имеет некоторые преимущества. Одним из важнейших преимуществ является вычисление параметров управления на основе физических факторов путем подбора подходящих значений параметров жесткости и демпфирования. Несмотря на тот факт, что данный подход учитывает магнитные подшипники как обычные механические подшипники и демпферы без учета остальных их многочисленных достоинств, было показано, что децентрализованное управление хорошо применимо к большому числу систем с АМП без значительных недостатков [2,6]. В большинстве таких случаев АМП системы с обратной связью демонстрируют приемлемое выполнение и надежные свойства, отвечающие требованиям стандарта ISO 14839 для магнитных подшипников [2,9].

9. Ортогонализация системы уравнений с приведением в бинормальные

координаты.

Поскольку в системе уравнений присутствует управление, она является неконсервативной, а матрицы Dx и K несимметричны. Это значит, что если мы хотим диагонализи-

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

Приведём систему (30) к форме Коши. Введём вектор z, состоящий из вектора перемещений и его производных:

г = {у q}Г

где вектор производных от вектора перемещений:

-си

Cистему (30) можно представить как:

^ л}л

В матричном виде система:

Здесь матрицы коэффициентов и вектор нагрузки:

Р„

0 М М Б,

, Оо

-М о

о к

ё =

(34)

Через введённые выше обозначения система (33) примет вид:

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

г1= К Чх}Г (36)

Данное уравнение поучено заменой матрицы коэффициентов основного уравнения транспонированными матрицами и изменением знака. С раскрытыми матрицами:

Это уравнение эквивалентно системе 2-ого порядка:

Предположим, что общее решение несимметричной системы уравнений 1-ого порядка имеет вид:

г = вх'С ■ ф (40)

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

(X Ро +Оо)ф = о (41)

Здесь ф - собственный вектор системы (35); X - соответствующее собственному вектору собственное число. Оно определяется действительной частью а и мнимой частью ю :

Х = а + гю (42)

Собственные векторы ф составляют матрицу собственных векторов

Ф = {ф ... ф8 } = [8 х 8]. Если представить аналогично выражению (40) решение исход-

ной системы (30), то получим и = {и ••• и} = [4х 8]. Это матрица собственных векторов для исходной системы. Поскольку значения корней характеристического уравнения X не зависят от формы уравнений задачи, мы приходим к выводу, что матрицы и и Ф связаны соотношением:

Хи и

ф

(43)

(44)

Аналогично полагая для сопряжённой системы (37):

^ = ецС • у

И подставляя в сопряжённую систему (37) будем иметь:

[ц-рТ= 0 (45)

Здесь у - собственный вектор системы (37), а ц - соответствующее ему собственное число. Собственные векторы у составляют матрицу собственных векторов Г = {у1 ... у8} = [8 х 8]. Если представить аналогично выражению (44) решение исходной сопряжённой системы второго порядка (39), то получим V = {^ ••• V } = [4 х 8]. Аналогично формам и и Ф , между матрицами форм V и Г существует связь:

"—П • V V

Корни ц отличаются от корней X только знаком, поскольку если ц заменить на —X, то получим то же характеристическое уравнение, что и для сопряжённой системы. Однако соответствующие формы ф и у будут разными:

Г:

(46)

Г:

"—п • V" "Х- V"

_ V _ _ V _

(47)

Найдём условие ортогональности, используя введённые формы. Представим уравнение для сопряжённой системы в виде:

уТ (ц-Ро —Qo) = 0 (48)

Подставим решение в виде гт = ектСт ■ фт и г1и = ец* Сп • уи в основную и сопряжённую систему соответственно:

|(Хт • Ро + Оо ) фт = 0

|тТ "(ци • Ро — Оо ) = 0

Умножим первое уравнение слева на уТп, а второе справа на фт, сложив их и сократив слагаемые, получим:

(Хт +ци )• УПРоФт = 0 (50)

Поскольку X =—ци, получим:

(49)

ГТРофт = 0 при т * П Ы ОоФт = 0

УП (Хт • Ро + Оо ) Фт = 0 при т =

(51)

п

Последние соотношения аналогичны условиям ортогональности для систем с симметричными матрицами коэффициентов и называются условиями биортогональности. Матрица Г сопряжённой системы должна вместе с матрицей Ф приводить Р0 и ^ к диагональной форме. Теперь мы можем найти систему независимых между собой уравнений относительно бинормальных координат e = {^ ... г% }Т :

г = Фе (52)

Далее перейдём на индекс к . Подставляя данное соотношение в основную систему и умножая на утк слева, будем иметь:

Используя соотношение биортогональности, получим:

Вводя в рассмотрение обобщенный коэффициент С и обобщённую вынуждающую функцию О :

I Ск = тТ РФк = 2Хк ■ уТ МИк + уТ D1 Ик

1Ок = 1Тк § = уТк fU

получим конечную форму независимого уравнения для к-ой моды:

(55)

Располагая начальными условиями г (0) = г0 для системы в координатах q, становится возможным решение задачи Коши в биортогональных координатах с начальными условиями е (0) = е0 = Ф_1г0:

Где X - диагональная матрица, на диагонали которой стоят элементы вектора собствен-

ных частот Хк , а Сс - вектор из отношений

О,

к .

С Л

Из системы (57) получим решение в бинормальных координатах. Для перехода в физические координаты используем соотношение (52). Таким образом, получаем решение во временной области.

10. Учёт внешнего вязкого трения в системе

Введём матрицу коэффициентов вязкого трения Б2 для общего случая, то есть когда она не является симметричной:

Вязкое трение введём в систему, основываясь на данных работ [11,12,13]. Используя подход Коги (Caughey), матрицу диссипативных коэффициентов вязкого трения представим в виде суммы:

л ^

в2 = М^п (М_1Ку = пК +... + пк(М_1К)

г=1

(60)

Здесь п - коэффициенты пропорциональности. Индекс г изменяется от единицы до

числа степеней свободы системы для учёта четырёх резонансных пиков: г = 1,4. Найдём неизвестные п .

Приведём матрицу Б2 к диагональному виду. Поскольку для системы (30) условие

ортогональности нам не известно, а из соотношений биортогональности (51) невозможно выделить диагонализированную матрицу, то получим собственные формы для ортогона-лизации из консервативной части системы (30):

Здесь К8ут - симметричная часть матрицы жёсткостей:

Ку = 1 (K + Кт ) (62)

Решая задачу на собственные значения для системы (30), находим матрицу собственных векторов иа для консервативной системы. Тогда диагонализированная матрица коэффициентов вязкого трения:

02^ =

Перепишем правую часть (63) в виде с учётом (60):

и^и, = VI ( л1К+ I- «+к (м-1к;г) и.

Произведения и^ К (М_1К) ий представляют собой матрицы типа:

(63) (64)

Введём вектор диагонализированных коэффициентов вязкого трения иш , состоя-

щий из элементов на диагонали 02а. :

С учётом (65) и (66) выражение (63) принимает вид:

Отсюда находим вектор неизвестных коэффициентов демпфирования n :

n = Fd2diag

F = A 1

В индексной записи:

n =

Z Fd

2diag j

(69)

(70)

j=1

Компоненты вектора п имеют размерность соответственно с, с 3 , с 5 , с 7. Компоненты вектора , отвечающие углу поворота, измеряются в к г ■ м 2 ■ с _ Компоненты вектора

d2dia , отвечающие линейному смещению, измеряются в к г ■ с

2diag

j=1

Выражение для D2 с учётом найденного вектора n :

4 3 4 4 4

D2 = F^diagj + ... + K (М-1 К) Z Fjd2diag j = KZ(M K)'" ^ j j=1 j=1 i=1

Знаки суммирования поменяем местами:

4 -1 4 4 4

Z К (М1К)- Z Fijd2diag j =Z d2dag j Z ^K (M1) i=1 j=1 j=1 i=1

Введём следующее обозначение для внутренней суммы:

2 diag j

(71)

(72)

С учётом данного обозначения, окончательное выражение для матрицы Б :

D2 =Z R( j 1 d2dlag j = R(1)d2diag 1 + ... + R( 4 d

t( 41

2 diag 4

(74)

j=1

Матрица диссипативных коэффициентов при векторе производных координат будет определяться суммой Б и Б :

4

О = ^ + б2

(75)

Уравнения движения ротора (30) с учётом внешнего демпфирования (75) принимает

вид:

11. Приведение системы уравнений движения ротора к безразмерному

виду

Приведение системы к безразмерному виду позволит в дальнейшем, оперируя сравнительно небольшим количеством безразмерных параметров, получать наиболее общие результаты. Запишем в развёрнутом виде первые два уравнения системы (76):

Примем следующие величины в соответствующем масштабе:

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

Примем величины коэффициентов управления в следующем виде:

рА , р а к.

I,' I (79)

^А = ^ ^

к V

Приняты следующие масштабы: 5 - величина зазора, измеряется в м , т - масса ротора, измеряется в к г; Т - период вращения ротора, измеряется в с:

2л

Т=(80)

V - собственная парциальная частота консервативной системы при линейных колебаниях, измеряется в :

v =

к„

(81)

m

2 . „-1.

dcrl - критическое демпфирование по угловому смещению, измеряется в к г ■ м/ ■ с

dri = iJlÄky = г m2 • \к\12 = dr22 (82)

dcr2 - критическое демпфирование по линейному смещению, измеряется в к г ■ с _

dcr 2 = 2j m • |k| (83)

Перечислим безразмерные параметры:

в - безразмерное линейное смещение вдоль оси х (аналогично ^ - вдоль оси у). а - безразмерный параметр положения центра тяжести ротора (аналогично b, c, d -тильды опущены здесь и далее),

Ух, Уz - коэффициенты соотношения моментов инерции, х - безразмерное время,

p - безразмерная собственная частота консервативной системы при линейных колебаниях,

£. - безразмерные параметры демпфирования,

Ва - безразмерный коэффициент управления по скорости, Ра - безразмерный коэффициент управления по смещению. Разделив 1-ое уравнение системы на

2-ое на |к| '5 и введя безразмерные параметры, получим (значки тильды далее сняты):

1 т у /

Р+

+ ■

+ —

а' + Ъ1 к:Пл(ас + Ь^) + 21

а'+Ь1)

а1+Ь2)б1

1

Т

+ь-

)

о

т у у -/1 2/^у.т

+

М

+

д

й

+ Ь7)6

аЧЬ3

)

'лт^

1^1 к^) |5Д

р+

Е-

■а.+

+ ■

21

т\

N1

а'+Ъ1)

(

14 -=1 ^ Л14 ЧЗ

к,

I 14 ^ Л14 Ьч/

Л+

кр,

Тогда, учитывая комплекс соотношений (14) для частот, периода и токов управления для РБ-регулирования, а также, что к < 0 :

+

+

+

+

А ■А

2л

2л

* +61) (с,1+Ь1У^ 4 ^ 4' 5! („Чг,1) («^Г ' ® (а1+Ь1)

У

„ ^ +2 ^ (ллк:1!^!^

Е-

■а+

25

V"'* /.Я* V , 1 1Д--4 -Л + Ду

\ ]4 --1 ^ л14 44 /

р+

- £ =

М*

п+

'4

2п

£ +

2л:

+

+

ч

Для удобства записи разложим матрицу Б2 по нечётным и чётным индексам:

В2 = б213)+Б22*4)

<85)

(86)

На основе уравнений, приведённого выше, в матричном виде имеем следующую структуру уравнений:

Щ + (С+Ос + в^ + Ч + (Кк ■+ кс | ц =

<К7)

Используя обозначения для матриц перед вектором перемещений и вектором производных перемещений (29), (75) и (86), окончательно имеем:

Система уравнений сохранила свою структуру, но составляющие необходимо переписать с учётом приведения к безразмерному виду:

Я = {Р в -а

(89)

М = \Р 12л,

У12

(а2 + Ь2) 0 1

О =

Гр 2

0 0

У г У х12

0 0

У х12

( а2 + Ь2) 0 0 1

(90)

Ос =1 — А. с 1 2л

2л Д а2 + Ь2) ( ас + Ьё) ( а + Ь )

0 0 10

0 0 0 0

-10 0 0

0 0 0 0

(91)

0

(а2 + Ь2) (а2 + Ь2) (с + ё) 2 0 0

( ас + Ьё) ( а + Ь )

Б

(1 3) — О /2

2/^Л/УГ

(а2 + Ь2)

• А

0

0

(1 3)

(а2 + Ь2) ( а2 + Ь2) 0 ( с + ё ) 2

(92)

-• А1 3)5 а12 5

( а2 + Ь2) 12 (а2 + Ь2)

А23)52 а23,3)5

• А

(1 3)

13

(а2 + Ь2)

А43)52

• 3)5

а14 5

• А

(1 3)

• 3)5

• А

,(1,3)

( а2 + Ь2) 31 ( а2 + Ь2) 32 (а2 + Ь2) 33 ( а2 + Ь2)

• 3)5

•Л 34

А13)5

а!13)52

А13)5

А 13)52

(93)

Б

(2,4) _

= 2^2

(а2 + Ь2) 5

(а2 + Ь2) 52

К

,(2,4)

А

,(2,4)

А

( 2,4)

А

( 2,4)

(94)

2

2

1

1

2

к

88

1

(а + Ь ) 0

0

( а + Ь ) ( а2 + Ь2) 2

0 0

00

00 1 (а + Ь )

(а2 + Ь2) (а + Ь) 2

(95)

Кс = Ра

( ас + Ьё) (а + Ь )

/и =

(а2 + Ь2) (а2 + Ь2) (с + ё) 2

0 0

0 0

0 0

0 ( ас + Ьё) (а + Ь )

(а2 + Ь2) (а2 + Ь2)

0 ( с + ё) 2

Л ^а Л

И5 К (а2 + Ь2' )52 5

(96)

(97)

или

/и = К /в ^0 /Л}Г (98)

12. Анализ устойчивости системы

Линейная неоднородная система дифференциальных уравнений устойчива тогда и только тогда, когда устойчива соответствующая однородная система. Поэтому рассмотрим однородную часть системы (88):

Матрица Кс (а значит и К ) несимметричны вследствие несовпадения областей действия магнитных сил со стороны АМП и датчика, регистрирующего смещение ротора. Ввиду несимметричности может произойти потеря устойчивости. Решение системы уравнений (99) будем искать в виде [8]:

q = %вх' (100)

где - собственный вектор; собственное значение X представим в комплексной форме:

Х = а + /ю (101)

где а - показатель роста, ю - собственная частота ротора.

После подстановки решения (101) в систему (100) получим:

(X2 • М + Х-Б + К)я0 = 0 (102)

Условие нетривиального решения данной системы:

ёе* (X2 • М + Х- Б + К) = 0 (103)

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

т

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

Для нахождения критического значения шс подставим выражение для собственного значения X в форме:

Хс = тс (104)

в уравнение (103):

ёа (-ш2 • М + тс- Б + К) = 0 (105)

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

Чтобы провести анализ устойчивости системы без учёта внешнего демпфирования, достаточно из уравнения (87) исключить слагаемые, отвечающие за диссипацию: Б(1,3) и

б22'4) .

13. Результаты расчёта

Исследование динамики ротора на АМП проводилось в программном комплексе МАТЬАВ©. Моделирование выполнено с возможностью наблюдать поведение системы при варьировании параметров демпфирования, частоты вращения ротора и жёсткостей магнитного подвеса.

Исходные данные приведены в таблице ниже:

Таблица 3. Исходные данные для исследования динамики

Параметр Символ Размерность Значение

Плотность материала ротора Р кг ■ м3 4540

Длина ротора 1 м 22-10"3

Диаметр ротора ¿1 м 15,6-10"3

Расстояния от центра масс до подшипников а Ь м —9 ■ 10~3 9■10"3

Расстояния от центра масс до сенсоров с а м -10-Ю"3 10-Ю"3

Зазор между ротором и АМП 8 м 10-4

Амплитуда синусоидального возмущения А1 н 10-Ю"3

Частота синусоидального возмущения а*, рад ■ с"1 112,5

Отрицательная позиционная жёсткость АМП К Н -м"1 -103

Токовая жёсткость АМП к, Н - А"1 0,01

Безразмерные параметры демпфирования по угловой степени свободы - 0,03

Безразмерные параметры демпфирования по линейной степени свободы Ь, ^4 - 0,03

Частота вращения ротора а рад ■ с"1 700

Значения коэффициентов управления РА, ИА устанавливаем соотношениями (26).

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

методу Коги (Са^Ьеу) (60). Решая задачу на собственные значения для основной (35) и сопряжённой (37) систем с учётом внешнего вязкого трения, получаем новый набор собственных частот и форм.

Решение системы уравнений в бинормальных координатах (57) во временной области найдено при начальных условиях:

г0 ={0 0 0 0 0,01 -0,05 -0,01 0,05}г (106)

Отметим, что составляющие вектора г0 являются безразмерными физическими величинами, то есть необходим перевод вектора начальных условий в бинормальные координаты. Возврат из бинормальных координат обеспечивается соотношением (52). Систему уравнений следует решать на интервале безразмерного времени т такого, чтобы система вышла на устойчивое движение. В данной работе т = 15.

Ниже проведено сравнение состояний системы без учёта демпфирования и с его включением.

Из двух таблиц ниже следует тенденция к понижению мнимых частей собственных частот при вводе демпфирования, что согласуется с [11].

Таблица 4. Собственные частоты без учёта демпфирования

Ь1 Ь 2 Ь 3 Ь 4

-127,34 + 4,79-1 -1,41 + 0,05 ■ 1 -0,63 + 1,10-1 -0,63 + 1,10-1

Таблица 5. Собственные частоты с учётом демпфирования

Ь1 Ь 2 Ь 3 Ь 4

-165,04 + 4,77 ■ 1 -1,08 + 0,03 ■ 1 -1,15 + 0,53-1 -1,15 + 0,53-1

Динамика свободных затухающих колебаний проиллюстрирована на двух рисунках

ниже.

Рисунок 6. Колебания системы без учёта демпфирования в физических координатах д

Рисунок 7. Колебания системы с учётом демпфирования в физических координатах д

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

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

Рисунок 8. Колебания системы с учётом демпфирования и синусоидального возмущения по оси х в

физических координатах д

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

г 1 494 1рпе _ 14__ ~ 83

ХА ~ 0,018 ~

(107)

Здесь гООЕ - время решения системы уравнений через встроенную функцию МЛТЬЛВ, измеряется в с,

- время вычисления аналитического решения, измеряется в с.

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

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

Рисунок 9. Колебания системы с учётом демпфирования и синусоидального возмущения по оси х для

бинормальных координат ег (т) . . . е4 (т)

ш

и-!Н

^

о.

о

о ^

ф -О

г о. о

а.2

0.15 0.1 0 05 0

-П 05 -0.1 -0.15 -0.2

Колебания в бинормапьных координатах е(5)... е(8)

и

5 ю

Безразмерное время, т

е5 е7

V \ и и

V \\ | Л :| Л 7 V

Г

V У

15

Рисунок 10. Колебания системы с учётом демпфирования и синусоидального возмущения по оси х для

бинормальных координат е5(т) ... е8(т)

14. Выводы

В ходе данной работы была построена модель жёсткого ротора на активном магнитном подвесе с учётом ПД-управления. Составлены уравнения движения ротора с учётом внешнего трения по схеме пропорционального демпфирования Коги (Са^Ьеу). Система уравнений приведена к безразмерному виду. Проведён анализ устойчивости системы, а также выведен алгоритм решения в бинормальных координатах для временной области. Анализ динамики с привлечением метода бинормальных координат позволяет исследовать несимметричные системы, что существенно расширяет возможности анализа линейных систем.

Учёт внешнего демпфирования является сложной инженерной задачей. В рамках данной работы был введён один из самых простых вариантов учёта внешнего трения. Он не учитывает возможности появления нестационарности в процессе обтекания объекта. В настоящее время появляются новые методики оценки влияния воздействия (и не только демпфирующего) внешней среды на обтекаемое твёрдое тело [1,14].

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

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

Список литературы

1. Шамолин М.В. Динамические системы с переменной диссипацией: подходы, методы, приложения // Фундаментальная и прикладная математика. 2008. Т. 14, № 3. С. 3-237.

2. Овсянникова Е.Е, Богданова Ю.В., Гуськов А.М. Исследование влияния потока крови на динамику ротора искусственного желудочка сердца (ИЖС) на активных магнитных подшипниках (АМП) // Наука и образование. МГТУ им. Н.Э. Баумана. Электрон. журнал. 2015. № 9. С. 298-317. DOI: 10.7463/0915.0811659

3. Крупнин А.Е., Банин Е.П., Гуськов А.М. Исследование поведения потока крови в области спрямителя искусственного желудочка сердца аксиального типа // 27 Международная Инновационно-ориентированная конференция молодых ученых и студентов (МИКМУС-2015): тез. докл. (Москва, 2-4 декабря 2015 г.). М.: ИМАШ РАН, 2015. С. 59.

4. Иткин Г.П. Механическая поддержка кровообращения: проблемы, решения и новые технологии // Вестник трансплантологии и искусственных органов. 2014. Т. 16, № 3. С. 76-84. DOI: 10.15825/1995-1191-2014-3 -76-84

5. Schweitzer G., Maslen E.H. Ch. 2: Principle of Active Magnetic Susoension; Ch. 7: Dynamics of the Rigid Rotors; Ch. 8: Control of the Rigid Rotor in ABMs // In book: Magnetic bearings. Theory, Design, and Application to Rotating Machinery. Springer Berlin Heidelberg, 2009. P. 31-53; P. 167-183; P. 191-197.

6. Журавлев Ю.Н. Активные магнитные подшипники: теория, расчет, применение. СПб.: Политехника, 2003. 206 с.

7. Schweitzer G. Active magnetic bearings - chances and limitations. Режим доступа: http://www.mcgs.ch/web-content/AMB-chances_and_limit.pdf (дата обращения 10.03.2016).

8. Genta G. Part 3: Dynamics of Rotating and Reciprocating Machinery // In book: Vibration Dynamics and Control. Springer US, 2009. P. 579-629.

9. ISO 14839-1:2002. Mechanical vibration - Vibration of rotating machinery equipped with active magnetic bearings - Part 1: Vocabulary. 2014. 30 p. // ISO Standarts: website. Режим доступа:

http://www.iso.org/iso/home/store/catalogue_tc/catalogue_detail.htm?csnumber=25726 (дата обращения 20.04.2016).

10. Халфман Р.Л. Динамика: пер. с англ. / под ред. И.А. Маркузона. М.: Наука, 1972. 568 с. [Halfman R.L. Dynamics. London: Addison Wesley Published Company, 1972.].

11. Бидерман В.Л. Теория механических колебаний. М.: Высшая школа, 1980. 480 с.

12. Gavin H.P. Classical damping, non-classical damping and complex modes. Department of Civil and Environmental Engineering Duke University. Fall, 2016.

13. Hale S. Damping in Ansys/LS-Dyna. Клуб пользователей ANSYS: сайт. Режим доступа: http://cae-club.ru/sites/default/files/users/files/1298/damping_dyna.pdf. (дата обращения 13.03.2016).

14. Пашков Е.Н., Мартюшев Н.В., Зиякаев Г.Р., Кузнецов И.В. Стационарное вращение неуравновешенного ротора, частично заполненного жидкостью при действии сил внешнего трения // Современные проблемы науки и образования. 2012. № 6. Режим доступа: http://science-education.ru/ru/article/view?id=7825 (дата обращения 01.07.2016).

Science ¿Education

of the Baumau MSTU

The Influence of External Damping upon the Dynamics of Rigid Rotor on Magnetic Bearings

S.V. Skoryukov1, A.M. Guskov1'*

Science and Education of the Bauman MSTU, 2016, no. 08, pp. 202-233.

DOI: 10.7463/0816.0844628

Received: 04.07.2016

Revised: 18.07.2016

© Bauman Moscow State Technical Unversity

gouskoY_am'Smail.m 1Bauman Moscow State Technical University, Moscow, Russia

Keywords: mechanical heart support system (MHSS), ventricular assist device (VAD), active

magnetic bearings (AMB), rotordynamics, stability of motion, biorthogonal coordinates

The paper conducts a linear analysis of the dynamics of rigid rotor installed in the active magnetic bearings (AMB). A current PD-AMP control is involved as the most simple and reliable algorithm.

The equations of motion of rigid rotor are obtained. Specific features of rheological blood are not taken into account and the blood is presented as a Newtonian fluid. Particular attention is paid to the inclusion of dissipative terms that takes into account the effect of blood on the dynamics of the impeller. External damping is included into the system by the method of Caughey (Caughey, 1965) to allow the required number of resonance peaks of the mechanical system. The included external damping expands a subject area of the paper by Y.V. Bogdanova, A.M. Guskov, E.E. Ovsyannikova «The research of influence of blood upon the dynamics of artificial ventricle rotor on active magnetic bearings».

The equations of motion are reduced to dimensionless form. The paper conducts analysis of the rotor motion stability. Using a method of asymmetric systems allowed us to reduce the equations of motion to a diagonal form and solve the system of the equations. To analyse dynamics were used the unrelated bi-orthogonal coordinates.

Further researches suppose introducing an integral term of the PD control (PID controller), a circulation force from the effects of blood flow added to the model, and conducting a nonlinear analysis of the rigid rotor dynamics.

This work may be useful to engineers majoring in rotor dynamics and specialists-designers of AMB control systems.

References

1. Shamolin M.V. Dynamical systems with variable dissipation: approaches, methods, and applications. Fundamentalnaya i Prikladnaya Matematika, 2008, vol. 14, no. 3, pp. 3-237. (English version of journal: Journal of Mathematical Sciences, 2009, vol. 162, no. 6, art. no. 741- DOI: 10.1007/s 10958-009-9657-y ).

2. Ovsiannikova E.E, Bogdanova Yu.V., Guskov A.M. The Research of Influence of Blood upon the Dynamics of Artificial Ventricle Rotor on Active Magnetic Bearings. Nauka i obrazovanie MGTU im. N.E. Baumana = Science and Education of the Bauman MSTU, 2015, no. 9, pp. 298-317. DOI: 10.7463/0915.0811659

3. Krupnin A.E., Banin E.P., Gus'kov A.M. Blood flow investigation in straightener domain ofaxial left ventricular assist device. 27 Mezhdunarodnaya Innovatsionno-orientirovannaya konferentsiya molodykh uchenykh i studentov (MIKMUS-2015): tez. dokl. [Abstracts of the 27 International Innovation Conference of Young Scientists and Students (IICYSS-2015)], Moscow, December 2-4, 2015. Moscow, Russian Academy of Sciences, Institute of Machines Science named by A.A. Blagonravov Publ., 2015, p. 59.

4. Itkin G.P. mechanical circulatory support: problems, solutions and new directions. Vestnik Transplantology i Iskusstvennykh Organov = Russian Journal of Transplantology and Artificial Organs, 2014, vol. 16, no. 3, pp. 76-84. DOI: 10.15825/1995-1191-2014-3-76-84

5. Schweitzer G., Maslen E.H. Ch. 2: Principle of Active Magnetic Susoension; Ch. 7: Dynamics of the Rigid Rotors; Ch. 8: Control of the Rigid Rotor in ABMs. In book: Magnetic bearings. Theory, Design, and Application to Rotating Machinery. Springer Berlin Heidelberg, 2009, pp. 31-53; pp. 167-183; pp. 191-197.

6. Zhuravlev Yu.N. Aktivnye magnitnye podshipniki: teoriya, raschet, primenenie [Active magnetic bearings: theory, calculation, application]. St. Petersburg, Politekhnika Publ., 2003. 206 p. (in Russian).

7. Schweitzer G. Active magnetic bearings - chances and limitations. Available at: http://www.mcgs.ch/web-content/AMB-chances_and_limit.pdf , accessed 10.03.2016.

8. Genta G. Part 3: Dynamics of Rotating and Reciprocating Machinery. In book: Vibration Dynamics and Control. Springer US, 2009, pp. 579-629.

9. ISO 14839-1:2002. Mechanical vibration - Vibration of rotating machinery equipped with active magnetic bearings - Part 1: Vocabulary. 2014. 30 p. ISO Standarts: website. Available at:

http://www.iso.ors/iso/home/store/catalogue tc/catalogue detail.htm?csnumber=25726 , accessed 20.04.2016.

10. Halfman R.L. Dynamics. London, Addison Wesley Published Company, 1972. (Russ. ed.: Halfman R.L. Dinamika. Moscow, Nauka Publ., 1972. 568 p.).

11. Biderman V.L. Teoriya mekhanicheskikh kolebaniy [Mechanical vibrations theory]. Moscow, Vysshaya shkola Publ., 1980. 480 p. (in Russian).

12. Gavin H.P. Classical damping, non-classical damping and complex modes. Department of Civil and Environmental Engineering Duke University. Fall, 2016.

13. Hale S. Damping in Ansys/LS-Dyna. ANSYS Club: website. Available at: http://cae-club.ru/sites/default/files/users/files/1298/damping_dyna.pdf , accessed 13.03.2016.

14. Pashkov E.N., Martyushev N.V., Ziyakaev G.R., Kuznetsov I.V. Stationary unbalanced rotor rotation which has been partially filled with liquid at external friction forces action.

Sovremennye problemy nauki i obrazovaniya = Modern problems of science and education, 2012, no. 6. Available at: http://science-education.ru/ru/article/view?id=7825 , accessed 01.07.2016. (in Russian).