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

УДК 62-251: 514.853

БОТ: 10.18698/2308-6033-2018-10-1815

Кинематически точное разделение большого поворота на осевой и поперечный в задачах роторной динамики

© Ф.Д. Сорокин, Х. Чжан МГТУ им. Н.Э. Баумана, Москва, 105005, Россия

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

Ключевые слова: большие повороты, вектор Эйлера, тензор поворота, тензор Жилина, композиция поворотов вращающихся элементов, роторная динамика

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

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

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

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

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

Углы Эйлера не применимы при нулевом угле нутации, углы

%

Крылова не позволяют перейти рубеж —, для вектора Родрига поворот ограничен величиной п, для вектора Эйлера — величиной 2п. Та-

%

ким образом, интервалы углов 0... %; 0... —, 0.. .2% для очень многих

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

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

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

п

может достигнуть величины —.

Описание больших поворотов ротора с использованием вектора Эйлера. Максимальным преимуществом и наибольшей наглядностью, по нашему мнению, обладает способ описания больших поворотов ротора с помощью вектора Эйлера [6-10]. Векторное описание поворотов ротора основано на теореме Эйлера о том, что произвольная комбинация пространственных поворотов ротора эквивалентна его одному плоскому повороту. Именно вектор Эйлера задает этот плоский поворот. Его направление совпадает с осью поворота, а длина равна углу поворота (рис. 1).

С вектором Эйлера связаны следующие тензорные функции векторного аргумента [6]:

L (S) = E cos 0 + SS + ^ Sx E;

е2

е

B (S) = E ^ + SS + Sx E,

е

е3

е2

где & — вектор Эйлера; е — длина вектора Эйлера (см. рис. 1); L(&) — функция, по которой вычисляют тензор поворота ротора по заданному вектору Эйлера; B(&) — функция, по которой вычисляют тензор Жилина по заданному вектору Эйлера; E — единичный тензор; && — диадное произведение векторов; — кососимметрич-ный тензор с сопутствующим вектором

Отметим, что понятие «тензор Жилина» не относится к общепринятым, например в работе [11] использована матрица с названием tangent operator, соответствующая тензору B . Хотя П.А. Жилин и не является автором этого тензора, именно после публикации его работ тензор B превратился в повседневный инструмент очень многих

исследователей. Поэтому использование фамилии Жилина в названии тензора В вполне оправданно.

Тензор поворота хорошо известен, его назначение понятно из названия, он «поворачивает» радиус-вектор.

Тензор Жилина связывает вектор угловой скорости © с производной вектора Эйлера

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

С помощью уравнения (2) можно вычислить криволинейный интеграл по годографу изменения вектора Эйлера.

Следует отметить, что приращение вектора Эйлера не является физическим поворотом некоторого узла, и это видно из уравнения (2). Фактически, момент должен совершать работу на произведении В (д)- <д, поэтому «физическим» малым поворотом 5ф в случае изменения вектора Эйлера на малый вектор Ад является произведение

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

Приложение векторно-тензорной теории больших поворотов ротора к построению различных конечных элементов [11-13] в значительной степени основано на соотношении (3).

Вектор Эйлера в сочетании с тензором поворота и тензором Жилина дает полноценное геометрически нелинейное описание вращения.

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

ш = В (д)- —

(1)

(2)

5ф = В (д)-Ад,

(3)

поворота ротора вокруг фиксированной продольной оси будут действовать обычные правила кинематики плоского движения из теоретической механики. Поэтому полный поворот ротора с вектором Эйлера д удобно представить как комбинацию двух последовательных поворотов: первый — вокруг оси ротора е на угол ф, второй — вокруг оси, перпендикулярной оси ротора, с вектором Эйлера у (рис. 2).

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

Рис. 2. Осевой (угол ф) и поперечный (вектор Эйлера у) повороты

L д = L у - L ф ;

Ь д= Ь (д); Ьу= Ь (у ); Ьф = Ь (фе).

(4)

Тензоры поворота дифференцируются согласно следующему правилу [6]:

= ш х Ь

<Ьф

= шфхЬф;

= шу хЬу,

(5)

где t — время; © — полная угловая скорость; ©ф — угловая скорость от вращения ротора вокруг продольной оси; ©у — угловая скорость, связанная с изменением вектора у.

Из формул (4) и (5), как показано в работе [6], следует правило комбинирования угловых скоростей при наложении двух последовательных поворотов:

ш = ш у + Ь у

ш

ф-

(6)

Таким образом, полный поворот удается точно определить с использованием разделенных векторов Эйлера фе и у: первый поворот ротора вокруг неподвижного орта е на угол ф, второй — вокруг поперечной оси, проходящей через вектор у на угол |у|.

Поскольку первый поворот — плоский, для него справедливы элементарные формулы кинематики плоского движения:

ш ф = ©фе;

< ф

(7)

©ф =

ф Л

Для второго поворота, согласно уравнению (1), угловую скорость находят с помощью тензора Жилина:

ю 7 = В • ^;

7 7 Ж (8)

В 7= В (Г )•

Подстановка уравнения (6) в выражения (8) с учетом выражений (7) позволяет выразить производную вектора у через угловые скорости:

^ = < (ю-шфЬу.е). (9)

Проецирование соотношения (9) на продольную ось ротора с учетом условия ортогональности векторов е и 7 приводит к тождеству

е • В-1 (ю-ШфЬу-е) = 0. (10)

Из тождества (10) следует выражение юф через полную угловую скорость:

е • В-1 • ю

Шф=е^. (11)

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

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

Рис. 3. Положение ротора в начальный момент времени

Блок пружин состоит из трех пружин растяжения — сжатия с жесткостью С1 и двух пружин кручения с жесткостью С2. Ротор представляет собой тонкий диск массой т и невесомый абсолютно жесткий стержень длиной I. В начальный момент времени, когда

ротор расположен горизонтально, ему задается начальная угловая скорость ©0 вокруг продольной оси.

Численный расчет был проведен для следующих исходных данных:

Диаметр диска Д м ................................................................................................................................................................0,5

Масса ротора т, кг..................................................................................................................................................................5

Расстояние от центра масс до центра

блока пружин I, м ..........................................................................................................................................................0,75

Ускорение свободного падения g, м/с2..................................................................................................9,81

Коэффициенты жесткости пружин С1, Н/м ..................................................................................2000

С2, Нм/рад ................................................................................................................................................................................500

Начальная угловая скорость ю0, рад/с ....................................................................................................50

Рассмотрим силы и моменты, действующие на ротор в процессе движения (рис. 4). Блок пружин создает в точке а силу F, направленную противоположно радиусу-вектору точки а, и момент М, направленный противоположно вектору у. Кроме реакций пружин на ротор действует сила тяжести mg, приложенная в центре тяжести диска (точка с).

Рис. 4. Силы и моменты, действующие на ротор в процессе движения

Динамику ротора описывает следующая система дифференциальных уравнений:

drc _ ~сй ~ Ус; d ф

и=иф ;

dy

dt

сЫ^ F

— = В"1 -(о-ШфЬ у- е); F

_-+ё;

т

_(а - Гс )Х F + М,

(12)

dt т dK

dt

где гс — радиус-вектор центра масс диска; ус — вектор скорости центра масс; К — вектор кинетического момента; га — радиус-вектор левого торца ротора.

При составлении системы (12) были использованы кинематические соотношения (7), (9), (11). Остальные уравнения следуют из обычных уравнений кинематики и динамики твердого тела, рассматриваемых в теоретической механике.

Система (12) замыкается алгебраическими уравнениями, связывающими вектор кинетического момента с полной угловой скоростью, представлением радиуса-вектора га через радиус-вектор гс и тензор поворота Ьу, а также соотношениями упругости:

—1 т

Ю = Jс • Кс; Iс = • Iс0 • ьа; Га = Гс — L у • е1;

р = — сЛ;

м = —С2 у,

где Лс — тензор инерции вращения в текущий момент времени; Лс0 — тензор инерции вращения в начальный момент времени

Гт то2 (

Лс 0 = (Е + ее) .

I 16 )

Система (12) была проинтегрирована со следующими начальными условиями:

Гс (0) = /е; ф(0) = 0; у(0) = 0; V с (0) = 0; К с (0) = I с0 • еш0.

Для численного интегрирования применялся традиционный метод Рунге-Кутты 4-го порядка с автоматическим выбором шага [14]. По результатам расчета построена траектория движения центра масс (рис. 5).

Установлена зависимость осевого поворота ф от времени, сопоставимая с классическим линейным законом (рис. 6). Из графика (см. рис. 6) видно, что зависимость ф(^) является нелинейной, хотя можно наблюдать весьма малое отклонение от линейного закона.

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

0,10

-0,16

-0,08 -0,04

0,4 г, м

Рис. 5. Траектория центра масс диска: а — пространственная; б — в плоскости уог

ф, рад

Рис. 6. Зависимость осевого поворота ф от времени (а) и его отклонение (б)

от линейного закона

Уу, рад 0,08

0,04 0

-0,04 -0,08

0

у2>рад

д лЛ А 1к

V | 1 г V/ | 1

2 з

а

с

-0,025 -0,050 -0,075 -0,100 -0,125 -0,150

0

2 3 4 г, с б

Рис. 7. Зависимость проекций вектора Эйлера у от времени на ось у (а), на ось г (б)

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

СОф-Юо> рад/с 0,02

-0,02

-0,04

Рис. 8. Отклонение осевой угловой скорости от начальной

Из данных на рис. 8 следует, что осевая угловая скорость изменяется во времени, но это отклонение весьма мало (наибольшее относи-

Ю™ - Ю0

не превышает 0,1 %).

тельное отклонение

Юг

Проверка результатов по закону сохранения энергии. Контроль численного решения осуществлялся проверкой выполнения закона сохранения энергии.

Кинетическая энергия в текущий и начальный моменты времени равна:

г 1 ^ 1 т

1 =— ш\г - ус +—ю• Тс • ю,

2 с с 2 с

Тп = Т(0) = |юо • Т

с0 - Ю0.

Потенциальная энергия состоит из энергии деформаций пружин и энергии диска в поле сил тяжести:

и = 1СЛ2 +1С2У 2 - - Гс

Установлена зависимость энергии от времени, причем в различных сочетаниях (рис. 9).

Полная энергия сохраняется с высокой точностью (тах |Т - Тп + и| =

= 2 -10-6 Дж), которая подтверждает отсутствие не только ошибок в представленных уравнениях, но и критических значений при использовании предлагаемого приема разделения поворота на осевой и поперечный.

Энергия, Дж

-2,0

Рис. 9. Зависимость комбинаций кинетической, потенциальной и суммарной энергии от времени:

1 — Т - Т0; 2 — Т - Т0 + и; 3 — и

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

ЛИТЕРАТУРА

[1] Челомей В.Н., ред. Вибрации в технике. Том 3. Колебания машин, конструкций и их элементов. Москва, Машиностроение, 1980, 544 с.

[2] Бранец В.Н., Шмыглевский И.П. Введение в теорию бесплатформенных инерциальных навигационных систем. Москва, Наука, 1992, 280 с.

[3] Bremer H. Elastic multibody dynamics: a direct Ritz approach. Dordrecht, Springer, 2008, 464 р.

[4] Журавлев В.Ф. Основы теоретической механики. 2-е изд., перераб. Москва, Издательство Физико-математической литературы, 2001, 320 с.

[5] Попов В.В., Сорокин Ф.Д., Иванников В.В. Разработка конечного элемента гибкого стержня с раздельным хранением накопленных и дополнительных поворотов для моделирования больших перемещений элементов конструкций летательных аппаратов. Труды МАИ, 2017, № 92.

URL: http ://trudymai.ru/published. php?ID=76832

[6] Жилин П.А. Векторы и тензоры второго ранга в трехмерном пространстве. Санкт-Петербург, Нестор, 2001, 276 с.

[7] Жилин П.А. Рациональная механика сплошных сред. Санкт-Петербург, Издательство Политехн. ун-та, 2012, 584 с.

[8] Rankin C.C., Brogan F.A. An element independent corotational procedure for the treatment of large rotation. Journal of Pressure Vessel Technology-Transactions of the ASME, 1986, vol. 108 (2), pp. 165-174.

[9] Crisfield M.A. Nonlinear finite element analysis of solid and structures. 2nd edition. Chichester, John Wiley&Sons, 2012, 544 p.

[10] Елисеев В.В., Зиновьева Т.В. Механика тонкостенных конструкций. Теория стержней. Санкт-Петербург, Издательство Политехн. ун-та, 2008, 95 с.

[11] Geradin M., Cardona A. Flexible Multibody Dynamics. A Finite Element Approach. Chichester, John Wiley&Sons, 2001, 340 p.

[12] Felippa C.A. A Systematic Approach to the Element-Independent Corotational Dynamics of Finite Elements. Department of Aerospace Engineering Sciences and Center for Aerospace Structures. Boulder, University of Colorado, 2000, 42 p.

[13] Felippa C.A., Haugen B. A Unified Formulation of Small-Strain Corotational Finite Elements: I. Theory. Computer Methods in Applied Mechanics and Engineering, 2005, vol. 194, pp. 2285-2335.

[14] Дьяконов В.П. Компьютерная математика. Теория и практика. Москва, Нолидж, 2001, 1296 с.

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

Ссылку на эту статью просим оформлять следующим образом: Сорокин Ф.Д., Чжан Х. Кинематически точное разделение большого поворота на осевой и поперечный в задачах роторной динамики. Инженерный журнал: наука и инновации, 2018, вып. 10. http://dx.doi.org/10.18698/2308-6033-2018-10-1815

Сорокин Фёдор Дмитриевич - д-р техн. наук, профессор кафедры «Прикладная механика» МГТУ им. Н.Э. Баумана. Автор более 70 публикаций в отечественных и международных научных журналах. Область научных интересов: роторная динамика. e-mail: sorokin_fd@mail.ru

Чжан Хао — аспирант кафедры «Прикладная механика» МГТУ им. Н.Э. Баумана. e-mail: zhang274234111@yandex.ru

Kinematically accurate separation of the large rotation into axial and transverse...

Kinematically accurate separation of the large rotation into axial and transverse in problems of rotor dynamics

© F.D. Sorokin, H. Zhang Bauman Moscow State Technical University, Moscow, 105005, Russia

In order to overcome the problem of "singular points", a method has been developed for the kinematically accurate separation of a large rotation into an axial (scalar) and a transverse (Euler vector). The proposal is based on the fact that in the problems of the rotor dynamics of machines consisting of shafts, gears, bearings, etc., the transverse rotation never reaches a value of 2n (a critical value for the Euler vector). The axial rotation is not limited in any way. Differential equations of the dynamics of rotational motion are also divided into axial (scalar) and transverse (vector). The above example shows that the constructed system of differential equations can be easily integrated by standard numerical methods to very large total rotations without any restrictions, provided that the transverse rotation is limited to 2 л. The control of the results is checked by observing the law of conservation of total energy.

Keywords: large rotations, Euler vector, rotation tensor, Zhilin tensor, composition of rotations, rotor dynamics

REFERENCES

[1] Chelomei V.N., ed. Vibrations in technology. Volume 3. Oscillations of machines, structures and their elements. Moscow, Mashinostroenie Publ., 1980, 544 p.

[2] Branets V.N., Shmyglevsky I.P. Introduction to the theory of freeform inertial navigation systems. Moscow, Nauka Publ., 1992, 280 p.

[3] Bremer H. Elastic multibody dynamics: a direct Ritz approach. Dordrecht, Springer, 2008, 464 р.

[4] Zhuravlev V.F. Fundamentals of theoretical mechanics. 2nd ed., revised. Moscow, Publishing House of Physical and Mathematical Literature, 2001, 320 p.

[5] Popov V.V., Sorokin F.D., Ivannikov V.V. Trudy MAI — Proceedings of the MAI, 2017, no. 92. Available at: http://trudymai.ru/published.php?ID=76832

[6] Zhilin P.A. Vectors and second-rank tensors in three-dimensional space. St. Petersburg, Nestor Publ., 2001, 276 p.

[7] Zhilin P.A. Rational Mechanics of Continuous Media. St. Petersburg, Polytechnic University Publ., 2012, 584 p.

[8] Rankin C.C., Brogan F.A. An element independent corotational procedure for the treatment of large rotation. Journal of Pressure Vessel Technology-Transactions of the ASME, 1986, vol. 108 (2), pp. 165-174.

[9] Crisfield M.A. Nonlinear finite element analysis of solid and structures. 2nd edition. Chichester, John Wiley&Sons, 2012, 544 p.

[10] Eliseev V.V., Zinovieva T.V. Mechanics of thin-walled structures. Theory of rods. St. Petersburg, Politechnical Uniersity Publ., 2008, 95 p.

[11] Geradin M., Cardona A. Flexible Multibody Dynamics. A Finite Element Approach. Chichester, John Wiley&Sons, 2001, 340 p.

[12] Felippa C.A. A Systematic Approach to the Element-Independent Corotational Dynamics of Finite Elements. Department of Aerospace Engineering Sciences and Center for Aerospace Structures. USA, University of Colorado, 2000, 42 p.

[13] Felippa C.A., Haugen B. A Unified Formulation of Small-Strain Corotational Finite Elements: I. Theory. Computer Methods in Applied Mechanics and Engineering, 2005, vol. 194, pp. 2285-2335.

[14] Dyakonov V.P. Kompyuternaya matematika. Teoriya i praktika [Computer mathematics. Theory and practice]. Moscow, Noledge/Knowledge Publ., 2001, 1296 p.

Sorokin F.D., Dr. Sc. (Eng.), Professor, Department of Applied Mechanics. Bauman Moscow State Technical University. Author of over 70 publications in national and international journals. e-mail: sorokin_fd@mail.ru

Zhang Hao, Ph.D. of the Department of Applied Mechanics, Bauman Moscow State Technical University. e-mail: zhang274234111@yandex.ru