Известия РАН. Механика твердого тела, 2020, № 2, стр. 40-58

1. Введение. Большая часть проблем в инженерной практике, связанных с вибрацией вала или ротора, возникает вследствие разрушительного влияния на систему статора. Оригинальность данной работы в том, что для определения возникающих при контакте сил нормального давления и сухого/полусухого трения используются гармонические перемещения, скорости и ускорения, т.е. уравнение нормальных колебаний, эквивалентной системы “гибкий ротор – податливый статор”. Рассматривается как система “гибкий ротор – упругие подшипники – фундамент”, кратко называемая исходной роторной, так и система “гибкий ротор – упругие подшипники – жесткий статор – упругие опоры – фундамент”, которая просто называется исходной роторно-статорной. В эквивалентной системе непосредственный контакт между гибким ротором и податливым статором представляется в виде дополнительной линейно упругой связи (пружины) между ними жесткостью k_k (контактная жесткость). Причем наряду с динамическими параметрами эквивалентной системы, называемыми параметрами контактной системы, активно используются динамические характеристики исходной системы при отсутствии в ней явного контакта ротора со статором, которые, в свою очередь, называются характеристиками бесконтактной системы.

Прежде исследователи определяли контактные силы либо динамикой податливого ротора, опертого на жесткий статор [1–5], либо жесткостью и демпфированием податливого статора [6–13], нормальная реакция которого в месте контакта имеет линейную [6–8] или нелинейную [9–11] или обобщенную в виде гистерезисной петли характеристику. В последнем случае учитываются как деформация статора, так и ее направление [12, 13], а также пластический эффект [13] на пути “нагрузка–разгрузка”. В некоторых работах [9–12] под ”статорной податливостью” понимается, в том числе перемещение статора на упруго-демпферных опорах.

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

Недостатками моделей с заданной упруго-диссипативной характеристикой статора N = являются эмпирические или теоретические коэффициенты. Например, в широко применяемой модифицированной модели Герца с вязкой компонентой, с устраненным “прилипанием” и коэффициентом восстановления, не зависящим от скорости сближения тел $N = {{k}_{k}}{{{\varepsilon }}^{{v}}} + {{d}_{k}}{{{\varepsilon }}^{w}}{\dot {\varepsilon }}$, где ${v} = 3{\text{/}}2$, $w = 1{\text{/}}4$, параметр k_k зависит от материалов тел и формы их поверхностей в области контакта, параметр d_k определяет значение коэффициента восстановления. В подобных моделях возникает проблема достоверного получения контактно-динамических характеристик системы, в частности максимальных ускорений ротора. Отсюда следует уделять отдельное внимание анализу корректности задания упругой компоненты характеристики статора и добиваться приемлемой точности принятой модели с помощью подбора параметров вязкой компоненты. Некорректный выбор параметров компонент может привести к физически неверным результатам, например, к значению коэффициента восстановления больше единицы.

2. Моделирование системы “гибкий ротор–податливый статор”. Для исследования динамики ротора, контактирующего со статором, необходимо построить колебательную модель, имеющую две степени свободы или больше, т.е. содержащую, по меньшей мере, несколько инерционных, упругих и демпфирующих элементов.

Роторная система или ситуация жесткозакрепленного статора. Любую осесимметричную роторную систему выгодно (поскольку просто, в том числе с точки зрения выполнения достаточно точных расчетов) представлять в виде комбинации трех конечных элементов. Конечные элементы ротора, в том числе их обобщенные силы и перемещения, показаны на рис. 1, где a – жесткий диск, b – балочный элемент, c – линеаризованная “короткая” опора.

a) Жесткий диск с сосредоточенными массой m_i, поперечным I_i и полярным I_0i моментами инерции. Уравнение двумерного¹¹ движения такого диска, вращающегося с постоянной угловой скоростью ω, в неподвижной системе отсчета и матричной форме [14]:

b) Вращающийся балочный элемент длиной l, радиуса r, с распределенными инерционно-упругими параметрами и внутренним трением гистерезисного типа. Уравнение движения такого балочного элемента в принятой форме [15, 16]:

Любой из элементов матриц этого уравнения является обобщенной силой в узле i, возникающей от единичного(й) ускорения (скорости) или перемещения в этом или i + 1 узле.

В данном общем случае элемент представляет собой балку Тимошенко, т.е. наряду с классическими процессами упругости (подматрицы [k]) и инерции прецессии, включающими поступательные перемещения (подматрицы [m_t]) и гироскопическое действие (подматрицы [g]), учитываются инерция вращения (подматрицы [m_r], добавляемые к [m_t]) и деформация сдвига (s-члены). Если отношение длины элемента к его диаметру получается больше единицы и меньше трех, то допустимо (с точки зрения сохранения высокой точности расчетов) пренебрежение влиянием деформации сдвига (использование балки Релея). Если ${1 \mathord{\left/ {\vphantom {1 {(2r)}}} \right. \kern-0em} {(2r)}} > 3$, то возможно рассмотрение без инерции вращения (т.е. балки Бернулли–Эйлера).

Принято, что балочный элемент обладает диссипацией, которая не зависит от частоты. Напряжение опережает соответствующую деформацию на угол ${\theta }$, который определяется материалом. Петля гистерезиса в плоскости напряжение–деформация повторяет эллипс, площадь которого пропорциональна sinθ и характеризует величину рассеянной энергии, так что $\sin {\theta } = {{{{{\eta }}_{{\text{h}}}}} \mathord{\left/ {\vphantom {{{{{\eta }}_{{\text{h}}}}} {\sqrt {1 + {\eta }_{{\text{h}}}^{2}} }}} \right. \kern-0em} {\sqrt {1 + {\eta }_{{\text{h}}}^{2}} }}$, где η_h – коэффициент гистерезисных потерь в материале балки.

c) Линеаризованная “короткая” опора, характеризующаяся уравнением вида:

т.е. с изотропной поперечной жесткостью и нулевой угловой жесткостью, а также с вязким демпфированием, выраженным коэффициентом ${{\eta }_{{v}}}$, и гистерезисными потерями в материале(ах) подшипника, выраженными коэффициентом η_hb.

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

Как правило, заметную роль играют только радиальные коэффициенты:

но даже в этой классической модели подшипника используются не все из них.

Существенные побочные коэффициенты жесткости k_xy, k_yx (относительно главных k_xx, k_yy) и вязкого демпфирования c_xy, c_yx (относительно c_xx, c_yy) характерны для подшипника жидкостного трения. Вместе с тем следует учитывать, что обычно такой подшипник значительно неконсервативен, т.е. его матрица жесткости (демпфирования) значительно асимметрична [17, 18]. Разность побочных коэффициентов жесткости Δk_ij = k_xy – k_yx определяет направление действия на ротор “циркуляционной” результирующей, что в случае немалости Δk_ij по сравнению с главными коэффициентами жесткости приводит к циркуляции (прецессии) ротора, причем при k_xy < k_yx к обратной, при k_xy > k_yx к прямой. Последнее непосредственно связано с возникновением явления динамической неустойчивости ротора на подшипниках жидкостного трения (в виде его значительной прямой асинхронной прецессии).

Здесь следует обратить внимание, что подшипник с равными главными коэффициентами жесткости (k_xx = k_yy) и малыми по сравнению с ними побочными коэффициентами жесткости (k_xy ≈ k_yx ≈ 0) называется изотропным. Подшипник с равными главными, а также с немалыми и почти одинаковыми побочными коэффициентами жесткости (k_xy ≈ k_yx ≠ 0) называется ортотропным. Подшипник с неравными главными (k_xx ≠ k_yy) и близкими к нулю или друг к другу побочными (k_xy ≈ k_yx ≈ 0 или k_xy ≈ k_yx ≠ 0) коэффициентами жесткости называется анизотропным, а подшипник со значительными и несовпадающими побочными (k_xy ≠ k_yx ≠ 0) – неконсервативным.

Также важно отметить, что место непосредственного контакта ротора со статором формально аналогично “короткой” опоре, лишенной вязкого трения, поэтому:

где η_hk – коэффициент гистерезисных потерь в материалах контактных поверхностей ротора и статора.

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

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

Число степеней свободы в плоскости n равно 2p, где p – число узлов, а (p – 1) – число балочных элементов. Векторы порядка $n \cdot 1$ образуются из совершающихся в плоскости поступательных и угловых перемещений узлов или из действующих в плоскости на некоторые узлы компонент внешних сил, среди которых могут быть и контактные. Подматрицы размера $n \cdot n$ означают сборку из подматриц или из констант конечных элементов.

Подматрицы инерции [M], гироскопических моментов [G] и жесткости [K] балочных элементов являются ленточными вида:

где под каждым двузначным индексом понимается номер строки и номер столбца матриц [m], [g] и [k].

Инерционная [M_d] и гироскопическая [G_d] подматрицы дисков, а также подматрица жесткости подшипников [K_b] диагональные вида ${{[ \ddots ]}_{{n \cdot n}}}$, поэтому их объединение с подматрицами балочных элементов может быть выражено как:

Здесь под каждым двузначным индексом понимается номер строки и номер столбца объединенных матриц, которые определяются тем или иным номером узла, константа C_i равна m_i, константа C_{i + 1} равна I_i или I_0i, константа C_j равна k_bj. В качестве примера можно посмотреть на рис. 2, где схематизируется произвольная роторная система и ее p узлов.

Константа С_k, равная k_k, добавляется в матрицу жесткости роторной системы ${{[K]}_{{n \cdot n}}}$ = [K] + [K_b], когда k-е или/и другое сечение гибкого ротора контактирует с упругой поверхностью неподвижного статора. Получающаяся матрица жесткости контактной (эквивалентной) роторной системы далее обособляется тильдой ${{[\tilde {K}]}_{{n \cdot n}}} = [K] + [{{K}_{b}}]$ + [K_k], где [K_k] – κ-заполненная диагональная матрица, κ – количество мест контакта.

Так как система осесимметричная, допустима более компактная форма (2.1) за счет объединения матричных уравнений

где ${{[{\text{M}}]}_{{n \cdot n}}} = [M] + [{{M}_{d}}]$, $[{\text{G}}] = [G] + [{{G}_{d}}]$.

В результате, определяя перемещения в плоскости xz действительными, а перемещения в плоскости yz мнимыми, т.е. рассматривая векторы комплексных обобщенных координат ${{\{ q\} }_{{n \cdot 1}}} = \{ {{q}_{x}}\} + {\text{i}}\{ {{q}_{y}}\} $ и сил ${{{\text{\{ F\} }}}_{{n \cdot 1}}} = \{ {{F}_{x}}\} + {\text{i}}\{ {{F}_{y}}\} $, будет:

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

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

При моделировании практически применяемой осесимметричной статорной подсистемы, т.е. состоящей из сравнительно жесткого и одновременно равножесткого корпуса, а также из нежестких и вместе с тем одинаковой жесткости элементов крепления к фундаменту, подходяще выбрать ее дискретное представление. Другими словами, в виде сосредоточенных в p + 1 узле массы m_s и момента инерции относительно поперечных осей I_s с линейными связями, в том числе c упруго-демпферными опорами суммарной в плоскости жесткостью k_s = k₁ + k₂ и общим для той же плоскости коэффициентом вязкого демпфирования η_vsk_s = η_vsk₁ + η_vsk₂. Произвольная роторно-статорная система и ее p + 1 узлов продемонстрирована на рис. 3.

Пусть координатами центра инерции статора являются x_s, y_s, φ_ys, φ_xs, и его поступательные/угловые колебания как отдельного твердого тела на податливых опорах малы, тогда уравнения его движения в поперечных направлениях и вокруг центра инерции (точки p + 1) примут вид:

Очевидно, такой подход обеспечивает наименьшее число необходимых степеней свободы роторно-статорной системы. Системными векторами и матрицами будут:

Причем матрица ${{[{\text{K}}]}_{{n + 2 \cdot n + 2}}}$ отличается от матрицы ${{[{\tilde {K}}]}_{{n + 2 \cdot n + 2}}}$ тем, что в первой k_k = 0.

3. Реакция осесимметричной системы “гибкий ротор–податливый статор” на контакт. Уравнения либо вида (2.1) с порядком матриц $2n \cdot 2n$ или $2(n + 2) \cdot 2(n + 2)$, либо вида (2.2) с порядком матриц $n \cdot n$ или $(n + 2) \cdot (n + 2)$ могут быть непосредственно использованы для определения реакции прецессирующего ротора на контактные силы. Пусть правую часть (2.1) или (2.2) образуют только компоненты контактных сил, приложенные к k-й точке исходной системы (рис. 4) и вызывающие регулярную (круговую установившуюся) прецессию:

Математические преобразования будут компактнее, если исходить из (2.2):

где $\{ {{f}_{N}}\} = {{\left\{ {..N..} \right\}}^{{\text{т}}}}$, $\{ {{f}_{T}}\} = {{\left\{ {..Т..} \right\}}^{{\text{т}}}}$, индексы $n \cdot n$ или $n + 2 \cdot n + 2$ опущены.

Как видно, это стандартная конечно-элементная формулировка системы “ротор–статор”. Исходный ротор на n подшипниках представлен комбинацией трех конечных элементов (жесткий диск, балочный элемент и линеаризованная вязкоупругая опора) и сосредоточенно нагружен (в k-й точке) силами нормального давления N и сухого трения T = μN.

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

Как следствие, логичным является следующее предложение – рассматривать эквивалентную систему – ротор на n + 1 опорах в виде комбинации тех же конечных элементов, в которой дополнительная опора в k-м узле есть нежесткий шарнир (свободное опирание на податливый статор), при этом разгруженную здесь же на величину результирующей Q тех же самых сил N и T = μN, а также с фазовым сдвигом поперечных перемещений (эквивалентной системы относительно исходной) на угол, равный углу трения, т.е. удовлетворяющий тождеству μ = tanγ (рис. 4). Изображены перемещения k-го узла исходной (темн.) и i-го узла эквивалентной (светл.) системы “гибкий ротор–податливый статор”, контактные силы (нормального давления N и сухого трения T = μN), проекции этих сил, нагружающие исходную систему, и их результирующая (Q), разгружающая эквивалентную систему.

Конечно-элементная формулировка эквивалентной системы “ротор–статор” с возросшей поперечной жесткостью (с бóльшими коэффициентами матрицы жесткости по сравнению с исходной системой), а также с противоположно направленными контактными силами:

Так как движение предполагается круговым, стационарным и обратным, то Ψ = = ‒Ωt и решения имеют вид:

где A_i = ρ_i = const, A_{i + 1} = φ_i = const, ${{\tilde {A}}_{i}} = {{{\tilde {\rho }}}_{i}} = {\text{const}}$, ${{\tilde {A}}_{{i + 1}}} = {{{\tilde {\varphi }}}_{i}} = {\text{const}}$, т.е. все амплитуды действительные, положительные и постоянные.

Подстановка этих равенств в вышеприведенные уравнения дает:

Если пренебречь гироскопическим влиянием и трением, то:

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

Это уравнение линейное, поэтому применим метод суперпозиции мод, т.е. применимы разложение реальных перемещений ротора без опоры {A} и с опорой $\{ \tilde {A}\} $ на статор по собственным формам колебаний свободного от контакта и опертого ротора, а также переход к соответствующим модальным перемещениям:

где ${{\{ {\text{u}}\} }_{{v}}}$ и ${{a}_{{v}}}$ – вектор безразмерных (относительных) узловых перемещений и их размер (модальная амплитуда) для ${v}$-й собственной частоты/формы бесконтактной системы, {ũ}_w и ã_w – модальные перемещения узлов для w-й собственной частоты/формы контактной системы, s – число рассматриваемых мод.

После подстановки (3.4) в матричные уравнения (3.2), преобразования каждого из них к приведенному виду, т.е. умножения слева на [M]^–1, а также после левого умножения их на “резонирующую” или характерную форму колебаний $\{ {\text{u}}\} _{i}^{{\text{т}}}$ и $\{ {\tilde {u}}\} _{j}^{{\text{т}}}$ (что представляет собой частотно-избирательный отклик колебательной системы на периодическое внешнее воздействие, который проявляется ростом амплитуды характерных для данной системы колебаний при приближении частоты внешнего воздействия к одной из собственных частот) будет:

где E – единичная матрица, ${\mathbf{D}} = {{[{\text{M}}]}^{{ - 1}}}[{\text{K}}]$, ${\mathbf{\tilde {D}}} = {{[{\text{M}}]}^{{ - 1}}}[{\tilde {K}}]$.

Из условия ортогональности собственных форм колебаний (которое можно объяснить как равенство нулю сил инерции i-й (j-й) формы колебаний на перемещениях ${v}$-й (w-й) формы колебаний) следует:

поэтому остаются лишь члены при i = ${v}$ и j = w:

Делая замены ${{\{ {\text{u}}\} }_{i}}{{a}_{i}} = {{\{ A\} }_{i}}$ и ${{\{ {\tilde {u}}\} }_{j}}{{\tilde {a}}_{j}} = {{\{ \tilde {A}\} }_{j}}$, можно получить то, что возможна альтернативная запись фундаментального уравнения (3.3):

где ${{\{ A\} }_{i}} = \{ \begin{array}{*{20}{c}} {{{A}_{{(1)}}}}& \ldots &{{{A}_{{(2k - 1)}}}}& \ldots \end{array}\} _{i}^{{\text{т}}}$ – известный вектор реальных узловых перемещений, соответствующих i-й собственной форме бесконтактной системы {u}_i = {u₍₁₎ ... ${{u}_{{(2k - 1)}}}\,\, \ldots \} _{i}^{{\text{т}}}$.

Вследствие нормирования по матрице масс (приведения к единичной инерционной матрице):

где λ_i и ${{{\tilde {\lambda }}}_{j}}$ – i-я и j-я круговые собственные частоты бесконтактной и контактной систем (здесь по определению невращающихся).

В результате существуют отношения:

Приравнивая, получится что:

Умножая справа на {u}_i, это равенство можно свести к зависимости модальной амплитуды для j-й роторной формы при опирании на статор от амплитуды модальных перемещений по i-й собственной форме ротора, свободного от контакта:

где ${{\Delta }_{{ji}}} = \{ {\tilde {u}}\} _{j}^{{\text{т}}}{{\{ {\text{u}}\} }_{i}}$, что можно назвать коэффициентом суперпозиции мод (наложения без взаимного влияния друг на друга i-й характерной моды ротора до его прикосновения к статору на j-ю характерную моду ротора после касания им статора).

Амплитуда модальных перемещений по i-й форме свободного ротора a_i определяется в задаче бесконтактной динамики соответствующего ротора, поэтому, например, как известно, при его неуравновешенности без учета внешнего трения:

где ω – скорость вращения, ω_i – i-я критическая скорость, {d} – вектор дисбалансов (кг ∙ м).

Зная вектор и амплитуду относительных (модальных) перемещений при контакте ротора со статором, можно определить реальные амплитуды (амплитуды реальных перемещений) на каждой “контактной” моде:

Суммируя реальные амплитуды на всех рассмотренных модах, можно найти амплитуды реальных контактных перемещений ротора:

Так как выборка зазора между ротором и статором – результат динамического прогиба бесконтактной системы, то вектор ${{\{ A\} }_{i}} + {{\{ \tilde {A}\} }_{j}}$ – полные реальные перемещения узлов контактной системы. В сущности, суммарные перемещения ${{\{ A\} }_{i}} + {{\{ \tilde {A}\} }_{j}}$ определяют радиусы орбит кругового движения различных сечений ротора при его постоянном контакте со статором, которые показательно соотносить с величиной радиального зазора δ = R – r, где r и R – радиусы поверхностей ротора и статора в месте их контакта.

Так как {f_N} > 0, то на основании (3.2) должно быть:

Замены (3.4), умножения слева на [M]^–1, а затем на $\{ {\text{u}}\} _{i}^{{\text{т}}}$, $\{ {\tilde {u}}\} _{j}^{{\text{т}}}$, вследствие (3.5), дают:

В результате (3.6) и объединения неравенств:

Отношение (3.7) описывает реакцию невращающейся недемпфированной системы “гибкий ротор–податливый статор” на установившийся контакт, что возможно при выполнении условия существования или устойчивости безотрывного движения ротора по статору (3.8).

В случае зависимости частот и форм собственных колебаний λ_ωi, {u_ω}_i и ${{{\tilde {\lambda }}}_{{{\omega }j}}}$, {ũ_ω}_j от скорости вращения (когда учитываются эффекты вращения: гироскопическое действие, внутреннее, внешнее и контактное трение) “гироскопический вклад” и “вклад трения” могут быть проанализированы по отдельности. Это допустимо ввиду пренебрежимо слабого влияния трения на собственные частоты и формы колебаний.

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

Очевидно и ожидаемо, сила нормального давления {f_N} содействует силам упругости [K]{A} исходного ротора, который условно не опирается на статор, а силы трения (η_h[K] + η_hb[K_b] + Ω${{\eta }_{{v}}}$[K_b]){A} противодействуют силе (уравновешивают силу) контактного трения $\{ {{f}_{T}}\} $. Однако результирующая контактных сил {f_Q} = = =  противоположна восстанавливающим силам $[{\tilde {K}}]\{ \tilde {A}\} $ опирающегося (эквивалентного) ротора, т.е. сила Q по отношению к нему является возмущающей.

Если разложить бесконтактные {A} и контактные $\{ \tilde {A}\} $ реальные перемещения по собственным формам и перейти к модальным перемещениям, т.е. если подставить (3.4) в (3.9), а затем умножить полученные уравнения слева на [M]^–1, а также на характерные формы $\{ {\text{u}}\} _{i}^{{\text{т}}}$ и $\{ {\tilde {u}}\} _{j}^{{\text{т}}}$, то:

где S = [M]^–1[K], S_b = ${\mathbf{S}}_{b}^{*}$ = [M]^–1[K_b], S_k = [M]^–1[K_k], при колебаниях статора как твердого тела ${\mathbf{S}}_{b}^{*} \ne {{{\mathbf{S}}}_{b}}$, так как соответствующая матрица вязкого сопротивления (2.3) изменяется.

Собственные векторы динамических матриц ортогональны, поэтому выполняются равенства (3.5), (3.6) и:

где ${{\tilde {s}}_{j}}$, ${{\tilde {s}}_{{bj}}}$ или $\tilde {s}_{{bj}}^{*} = \{ {\text{u}}\} _{i}^{{\text{т}}}{\mathbf{S}}_{b}^{*}{{\{ {\text{u}}\} }_{i}}$, ${{\tilde {s}}_{{kj}}}$ (s_i, s_bi или $s_{{bi}}^{*}{\text{ = }}\{ {\tilde {u}}\} _{j}^{{\text{т}}}{\mathbf{S}}_{b}^{*}{{\{ {\tilde {u}}\} }_{j}}$) – собственное число вращающихся деталей (ротора), подшипников или опор в том числе статорных, и контакта соответственно для  j-й (i-й) формы колебаний контактной (бесконтактной) системы.

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

где ${{{\varepsilon }}_{{\eta }}} = {{{\eta }}_{{\text{h}}}}{{s}_{i}} + {{{\eta }}_{{{\text{h}}b}}}{{s}_{{bi}}} + {{{\eta }}_{{v}}}s_{{bi}}^{*}\Omega $, ${{{\tilde {\varepsilon }}}_{{\eta }}} = {{{\eta }}_{{\text{h}}}}{{\tilde {s}}_{j}} + {{{\eta }}_{{{\text{h}}b}}}{{\tilde {s}}_{{bj}}} + {{{\eta }}_{{{\text{h}}k}}}{{\tilde {s}}_{{kj}}} + {{{\tilde {\eta }}}_{{v}}}\tilde {s}_{{bj}}^{*}\Omega $.

Исключение $\{ {{f}_{N}}\} $ и $\{ {{f}_{T}}\} $ дает:

Умножение справа на {u}_i приводит к:

Неизвестные амплитудные множители ${{\tilde {a}}_{j}}\cos {\gamma }$ и ${{\tilde {a}}_{j}}\sin {\gamma }$ найдутся из отношений:

Отсюда, вследствие малости произведения ${{{\varepsilon }}_{{\eta }}} \cdot {{{\tilde {\varepsilon }}}_{{\eta }}}$, а также квадратов ${\varepsilon }_{{\eta }}^{{\text{2}}}$ и ${\tilde {\varepsilon }}_{{\eta }}^{{\text{2}}}$, во-первых:

во-вторых:

т.е. неизвестная Ω определяется уравнением четвертой степени: где ${{{\varepsilon }}_{{\text{h}}}} = {{{\eta }}_{{\text{h}}}}{{s}_{i}} + {{{\eta }}_{{{\text{h}}b}}}{{s}_{{bi}}}$, ${{{\tilde {\varepsilon }}}_{{\text{h}}}} = {{{\eta }}_{{\text{h}}}}{{\tilde {s}}_{j}} + {{{\eta }}_{{{\text{h}}b}}}{{\tilde {s}}_{{bj}}} + {{{\eta }}_{{{\text{h}}k}}}{{\tilde {s}}_{{kj}}}$, ${{{\varepsilon }}_{{v}}} = {{{\eta }}_{{v}}}s_{{bi}}^{*}$, ${{{\tilde {\varepsilon }}}_{{v}}} = {{{\eta }}_{{v}}}\tilde {s}_{{bj}}^{*}$ (в случае неподвижного статора ${{s}_{{bi}}} = s_{{bi}}^{*}$ и ${{\tilde {s}}_{{bj}}} = \tilde {s}_{{bj}}^{*}$).

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

При отсутствии вязкого и гистерезисного демпфирования будет:

С трением будет:

где λ – та или иная круговая собственная частота бесконтактной или контактной недемпфированной/невращающейся системы, Δ – некоторая поправка для λ от трения.

Представление как |Ω| = λ(1 + ε), где ε = Δ/λ – малая поправка, и подстановка этого и приближенных выражений

в уравнение четвертой степени раскрывают введенную ε-поправку:

Как правило, гистерезисные потери в материалах ротора, подшипников и статора незначительно разнятся ${{{\eta }}_{{\text{h}}}} \approx {{{\eta }}_{{{\text{h}}b}}} \approx {{{\eta }}_{{{\text{h}}k}}} \approx {\eta }$, где η – коэффициент гистерезисных потерь в материалах системы.

В результате:

Как следствие:

При ${{\eta }_{{v}}}$ = 0:

Имеющаяся неоднозначность ${\lambda } = {{{\lambda }}_{i}}$ или ${\lambda } = {{{\tilde {\lambda }}}_{j}}$ устранима, если проанализировать аналитические выражения:

С приближением η к μ значения скорости обкатки Ω по модулю в первом случае неограниченно возрастают, во втором – уменьшаются без разрыва. Физически увеличение трения не может давать увеличение скорости, отсюда λ не может быть равной λ_i.

Таким образом, ε-добавка к ${{{\tilde {\lambda }}}_{j}}$, определяющая Ω, т.е. угловую скорость обратного безотрывного движения ротора по статору, вследствие трения:

или

Поскольку $\tilde {s}_{{bj}}^{*} = \{ {\tilde {u}}\} _{j}^{{\text{т}}}{{[{\text{M}}]}^{{--{\text{1}}}}}[{{K}_{b}}]{{\{ {\tilde {u}}\} }_{j}}$, то как при малом вязком демпфировании, так и при мягких опорах ротора и/или статора члены с множителем η_v будут близкими к нулю:

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

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

4. Заключение. В отличие от существующих представлений при построении модели “ротор–статор” для определения возникающих при контакте сил нормального давления и сухого/полусухого трения используются не коэффициенты контактной жесткости и контактного демпфирования, а вибрационные перемещения, скорости и ускорения эквивалентной системы, т.е. уравнение ее нормальных колебаний. Тем самым априори неизвестные контактные силы исходной системы исключаются из рассмотрения. Причем под исходной системой понимается n-опорный ротор, загруженный на величину результирующей Q контактных сил, а под эквивалентной – ротор на n+1 опорах с возросшей поперечной жесткостью (с бóльшими по сравнению с n-опорным коэффициентами матрицы жесткости), и разгруженный на такую же величину Q.

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

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

где |Ω| – абсолютная величина угловой скорости обратной прецессии ротора, ω – скорость собственного вращения ротора, r – роторный радиус и δ – зазор в месте контакта, λ_i – i-я круговая частота собственных колебаний свободного от контакта ротора, ${{{\tilde {\lambda }}}_{j}}$ – j-я круговая частота собственных колебаний ротора, опирающегося на статор.

Во-вторых, угловая скорость прецессии ротора при обкатке статора с проскальзыванием между ними близка к ${{{\tilde {\lambda }}}_{j}}$, а именно:

где η – коэффициент гистерезисных потерь в материалах системы, причем πη равно логарифмическому декременту затухания в системе, μ – коэффициент трения скольжения.

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

Работа выполнена при финансовой поддержке Минобрнауки России в рамках Федеральной целевой программы “Исследования и разработки по приоритетным направлениям развития научно-технологического комплекса России на 2014–2020 годы”, Соглашение о предоставлении субсидии № 14.607.21.0191 от 26.09.2017 г., проект RFMEFI60717X0191.