Features of the dynamics of a rotating shaft with nonlinear models of internal damping and elasticity

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

Влияние внутреннего демпфирования, по всей вероятности, впервые было достаточно полно описано в работе Кимпбалла [5] и продемонстрировано в экспериментах Ньюкирка [6]. В последующем вопросы устойчивости вращаю­щихся роторов при наличии линейного внутреннего демпфирования были рассмотрены в работе [7]. Причем основным эффектом, связанным с наличием внутреннего демпфирования, является самовозбуждение изгибных колебаний. В последние годы усовершенствуется подход к моделированию эффектов внутреннего демпфирования ­[8–11]. Тем не менее особенности поведения в закритической области после бифуркации ранее подробно не рассматривались. Ограничение поперечных колебаний вала возможно как за счет внешних устройств [12], так и при учете нелинейностей в законе упругого деформирования вала.

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

2. Расчетная схема. Рассматривается гибкий вал круглого постоянного поперечного сечения, вращающийся вокруг своей продольной оси с постоянной угловой скоростью $\omega$. Оба конца вала установлены в жестких (недеформируемых) опорах, обеспечивающих свободное вращение вала, но исключающих смещение и поворот его концевых сечений относительно поперечных осей. Вал имеет погонную массу $m_R$, которая равномерно распределена по его длине $l$. Гироскопическими эффектами пренебрегаем. Вал моделируется стержнем Бернулли–Эйлера с включением дополнительных членов, учитывающих нелинейно-вязкое трение в модели Кельвина–Фойхта и кубической нелинейностью в законе упругости.

Учитываются силы внешнего трения вала, пропорциональные его абсолютной поперечной скорости колебаний. Поперечные колебания вала будем рассматривать относительно неподвижной системы координат $Oxyz$ с началом на левой опоре (рис. 1). Орты $\overrightarrow{i},\overrightarrow{j},\overrightarrow{k}$ связаны с осями $Oxyz$.

Рис. 1. Расчетная схема вращающегося вала: 1 – траектория прецессии, 2 – направление вращения.

Рис. 2. Диаграмма Аргана в диапазоне скоростей вращения $\Omega \in \left[30;32\right]$.

3. Учет внутреннего демпфирования при распределенной массе вращающегося вала. В линейной теории для описания внутреннего демпфирования в случае одноосного напряженного состояния обычно используется модели Кельвина–Фойхта для упруго-вязких тел [13]:

$\sigma =E\left(\epsilon +T_V\dot{\epsilon }\right)$,

где $E$ – модуль упругости материала, $\epsilon$ – деформация, $T_V$ – характерное время (произведение $E{T}_V$ представляет собой коэффициент вязкости).

Для учета внутреннего демпфирования в материале вала с равномерно распределенной массой необходимо во вращающейся системе координат определить осевую деформацию в выбранной точке материальной среды. Обозначим кривизну вала $\overrightarrow{\kappa }\left(z,t\right)$. В неподвижной системе координат $Oxyz$, $z\in \left[\mathrm{0,}l\right]$, кривизна имеет вид $\overrightarrow{\kappa }\left(z,t\right)={\kappa }_s\left(z;t\right){\overrightarrow{i}}_s$, $s=\mathrm{1,}2$ и во вращающейся вместе с валом системе координат $O\tilde{x}\tilde{y}z$, $z\in \left[\mathrm{0,}l\right]$, кривизна $\overrightarrow{\kappa }\left(z,t\right)={\tilde{\kappa }}_s\left(z;t\right){\overrightarrow{\tilde{i}}}_s$, $s=\mathrm{1,}2$. Базисы ${\overrightarrow{i}}_s$ и ${\overrightarrow{\tilde{i}}}_s$, $s=\mathrm{1,}2$ – ортонормированы, т.е. ${\overrightarrow{i}}_s\cdot {\overrightarrow{i}}_t={\tilde{\overrightarrow{i}}}_s\cdot {\tilde{\overrightarrow{i}}}_t={\delta }_{st}$. В силу инвариантности вектора кривизны верно равенство¹:

 $\overrightarrow{\kappa }\left(z;t\right)={\kappa }_s\left(z;t\right){\overrightarrow{i}}_s={\tilde{\kappa }}_s\left(z;t\right){\tilde{\overrightarrow{i}}}_s\left(t\right);s=\mathrm{1,}2$              (3.1)

Заметим, что наряду с инвариантной записью необходимо различать мат­ричные отображения ${\left.\overrightarrow{\kappa }\right|}_{O\tilde{x}\tilde{y}z}=\tilde{\kappa }$ и ${\left.\overrightarrow{\kappa }\right|}_{Oxyz}=\kappa$, которые привязаны к различным базисам. Связь между координатами в различных базисах определяется ортогональной матрицей поворота $S\left(\phi \right)$, где $\phi =\omega t$ – угол поворота:

$\kappa =S\left(\phi \right)\tilde{\kappa }:S\left(\phi \right)=\left[\begin{array}{cc}\cos \left(\phi \right)& \sin \left(\phi \right)\\ -\sin \left(\phi \right)& \cos \left(\phi \right)\end{array}\right]$.

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

$\begin{array}{l}{\epsilon }_z\left(\tilde{x},\tilde{y},z;t\right){\overrightarrow{i}}_3=-\overrightarrow{p}\times \overrightarrow{\kappa }=\left(-\tilde{x}{\tilde{\kappa }}_2\left(z;t\right)+\tilde{y}{\tilde{\kappa }}_1\left(z;t\right)\right){\overrightarrow{i}}_3;\\ \overrightarrow{p}=\tilde{x}{\tilde{\overrightarrow{i}}}_1+\tilde{y}{\tilde{\overrightarrow{i}}}_2,\overrightarrow{\kappa }={\tilde{\kappa }}_1\left(z;t\right){\tilde{\overrightarrow{i}}}_1+{\tilde{\kappa }}_2\left(z;t\right){\tilde{\overrightarrow{i}}}_2.\end{array}$

Или в матричной форме во вращающейся системе координат:

$\begin{array}{l}{\epsilon }_z\left(\tilde{x},\tilde{y},z;t\right)=-\tilde{x}{\tilde{\kappa }}_2\left(z;t\right)+\tilde{y}{\tilde{\kappa }}_1\left(z;t\right)={\tilde{p}}^{T}R\tilde{\kappa }\text{;}\\ \tilde{p}=\left\{\begin{array}{c}\tilde{x}\\ \tilde{y}\end{array}\right\},R=\left[\begin{array}{cc}0& -1\\ 1& 0\end{array}\right],\tilde{\kappa }=\left\{\begin{array}{c}{\tilde{\kappa }}_1\\ {\tilde{\kappa }}_2\end{array}\right\}.\end{array}$

Если для материала рассматривать линейный закон Гука и использовать модель стержня Бернулли–Эйлера, то изгибающий момент для сечения с двумя осями симметрии $I_{\tilde{x}}=I_{\tilde{y}}=I_R$ описывается следующим выражением [14]:

$\overrightarrow{M}=\tilde{\overrightarrow{M}}=E{I}_R\tilde{\overrightarrow{\kappa }}=E{I}_R\overrightarrow{\kappa }$.

Отметим, что координаты точки среды $\left(\tilde{x},\tilde{y}\right)$ во вращающейся вместе с валом системе отсчета при вращении не меняются, исходя из этого, скорость деформации ${\dot{\epsilon }}_z\left(\tilde{x},\tilde{y},z;t\right)$ для фиксированной точки среды $\left(\tilde{x},\tilde{y}\right)=\text{const}$ в мат­ричной форме имеет вид:

 ${\dot{\epsilon }}_z\left(\tilde{x},\tilde{y},z;t\right)=-\tilde{x}{\dot{\tilde{\kappa }}}_2\left(z;t\right)+\tilde{y}{\dot{\tilde{\kappa }}}_1\left(z;t\right)={\tilde{p}}^{T}R\dot{\tilde{\kappa }}$.             (3.2)

Вектор угловой скорости поворота $\overrightarrow{\omega }=\omega {\overrightarrow{i}}_3=\dot{\phi }{\overrightarrow{i}}_3$ сечения вала определяется скоростью поворота подвижного базиса относительно неподвижного базиса, тогда как ${\dot{\overrightarrow{\tilde{i}}}}_1=\overrightarrow{\omega }\times {\overrightarrow{\tilde{i}}}_1,{\dot{\overrightarrow{\tilde{i}}}}_2=\overrightarrow{\omega }\times {\overrightarrow{\tilde{i}}}_2$.

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

 $\begin{array}{c}\dot{\overrightarrow{\kappa }}=\frac{\partial }{\partial t}\left({\kappa }_s{\overrightarrow{i}}_s\right)=\frac{\partial }{\partial t}\left({\tilde{\kappa }}_s{\overrightarrow{\tilde{i}}}_s\right),{\dot{\overrightarrow{\tilde{i}}}}_s=\overrightarrow{\omega }\times {\tilde{\overrightarrow{i}}}_s;\\ \frac{\partial }{\partial t}\left({\kappa }_s{\overrightarrow{i}}_s\right)={\dot{\kappa }}_s{\overrightarrow{i}}_s,\frac{\partial }{\partial t}\left({\tilde{\kappa }}_s{\overrightarrow{\tilde{i}}}_s\right)={\dot{\tilde{\kappa }}}_s{\overrightarrow{\tilde{i}}}_s+{\tilde{\kappa }}_s{\dot{\overrightarrow{\tilde{i}}}}_s;\\ \dot{\overrightarrow{\kappa }}={\dot{\tilde{\kappa }}}_s{\overrightarrow{\tilde{i}}}_s+\overrightarrow{\omega }\times \left({\tilde{\kappa }}_s{\overrightarrow{\tilde{i}}}_s\right)={\dot{\tilde{\kappa }}}_s{\overrightarrow{\tilde{i}}}_s+\overrightarrow{\omega }\times \overrightarrow{\kappa };\\ {\dot{\tilde{\kappa }}}_s{\overrightarrow{\tilde{i}}}_s=\dot{\overrightarrow{\kappa }}-\overrightarrow{\omega }\times \overrightarrow{\kappa },\end{array}$                       (3.3)

где ${\dot{\tilde{\kappa }}}_s$ – производная компонент вектора кривизны, связанная с фиксированным материальным сечением вала $\left(z,t\right)$, которая входит в выражение скорости деформации.

Тогда уравнение состояния для вращающегося вала (одномерного объекта) от переменных $\left\{\overrightarrow{\kappa },\dot{\overrightarrow{\kappa }}\right\}$ в соответствии с линейной гипотезой внутреннего демпфирования в модели Кельвина–Фойхта будет уравнением изгибающего момента $\overrightarrow{M}\left(\overrightarrow{\kappa },\dot{\overrightarrow{\kappa }}\right)$, которое для вала с одинаковыми изгибными жесткостями относительно осей $O\tilde{x}$, $O\tilde{y}$ имеет вид:

 $\overrightarrow{M}\left(\overrightarrow{\kappa },\dot{\overrightarrow{\kappa }}\right)=E{I}_R\left(\overrightarrow{\kappa }+T_V{\dot{\tilde{\kappa }}}_s{\tilde{\overrightarrow{i}}}_s\right)=E{I}_R\left[\overrightarrow{\kappa }+T_V\left(\dot{\overrightarrow{\kappa }}-\overrightarrow{\omega }\times \overrightarrow{\kappa }\right)\right]$,             (3.4)

где $T_V$ – характерное время (время релаксации) [10], полагаем, что $T_V\ll 2\pi /\left|\overrightarrow{\omega }\right|$.

Для записи инвариантного уравнения в матричной форме в неподвижном базисе $Oxyz,\forall \overrightarrow{a}=a{\overrightarrow{i}}_3,\overrightarrow{b}=b_s{\overrightarrow{i}}_s$ используется матрица $R$:

 $\begin{array}{l}\overrightarrow{a}\times \overrightarrow{b}=a{\overrightarrow{i}}_3\times \left(b_1{\overrightarrow{i}}_1+b_2{\overrightarrow{i}}_2\right)=a\left(b_1{\overrightarrow{i}}_2-b_2{\overrightarrow{i}}_1\right);\\ R:a{\overrightarrow{i}}_3\times \overrightarrow{b}\iff aRb,\forall c\in R^2,c^{T}Rc=0;\\ R=\left[\begin{array}{cc}0& -1\\ 1& 0\end{array}\right],R^T=-R,R{R}^T=E,E=\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right].\end{array}$                 (3.5)

Тогда векторное уравнение в матричной форме в неподвижном базисе принимает следующий вид:

 $M=E{I}_R\left(\kappa +T_V\dot{\kappa }-T_V\omega R\kappa \right)$.            (3.6)

Последнее слагаемое в (3.6), зависящее от угловой скорости $\left(E{I}_R{T}_V\omega R\kappa \right)$, характерно для вращающихся деформируемых твердых тел. В данном случае именно это слагаемое в изгибающем моменте приводит к появлению циркуляционных сил и учитывает механизм передачи энергии вращения вала в изгиб вала, когда изгиб в одной плоскости вызывает изгиб в перпендикулярной плоскости, пропорционально скорости вращения вала $\omega$ и коэффициенту вязкости $E{T}_V$.

Учет дополнительных нелинейных членов в модели Кельвина–Фойхта можно провести в форме Коссера через скорость вектора кривизны в системе координат ${\dot{\tilde{\kappa }}}_s{\tilde{\overrightarrow{i}}}_s$, $s=\mathrm{1,}2$, вводя дополнительные слагаемые в законе состояния с вращающейся средой. Отметим, что скалярное произведение во вращаю­щейся системе ${\dot{\tilde{\kappa }}}_s{\dot{\tilde{\kappa }}}_s$ (суммирование по немому индексу $s=\mathrm{1,}2$ ) с учетом (3.3) вычисляется следующим образом в инвариантной форме:

${\dot{\tilde{\kappa }}}_s{\dot{\tilde{\kappa }}}_s=\left({\dot{\tilde{\kappa }}}_s{\tilde{\overrightarrow{i}}}_s\right)\cdot \left({\dot{\tilde{\kappa }}}_t{\tilde{\overrightarrow{i}}}_t\right)=\left(\dot{\overrightarrow{\kappa }}-\overrightarrow{\omega }\times \overrightarrow{\kappa }\right)\cdot \left(\dot{\overrightarrow{\kappa }}-\overrightarrow{\omega }\times \overrightarrow{\kappa }\right)$.

Например, по индукции можно учесть нелинейный, кубический член в линейной модели Кельвина–Фойхта (3.4), добавляя кубическое слагаемое ${\dot{\tilde{\kappa }}}_s{\dot{\tilde{\kappa }}}_s\left(\dot{\overrightarrow{\kappa }}-\overrightarrow{\omega }\times \overrightarrow{\kappa }\right)$:

 $\begin{array}{l}\overrightarrow{M}\left(\overrightarrow{\kappa },\dot{\overrightarrow{\kappa }}\right)=E{I}_R\left[\overrightarrow{\kappa }+\left(T_V+T_{VV}{\dot{\tilde{\kappa }}}_s{\dot{\tilde{\kappa }}}_s\right){\dot{\tilde{\kappa }}}_t{\tilde{\overrightarrow{i}}}_t\right]=\\ =E{I}_R\left[\overrightarrow{\kappa }+\left(T_V+T_{VV}{\dot{\tilde{\kappa }}}_s{\dot{\tilde{\kappa }}}_s\right)\left(\dot{\overrightarrow{\kappa }}-\overrightarrow{\omega }\times \overrightarrow{\kappa }\right)\right]=\\ =E{I}_R\left[\overrightarrow{\kappa }+\left[T_V+T_{VV}\left(\dot{\overrightarrow{\kappa }}-\overrightarrow{\omega }\times \overrightarrow{\kappa }\right)\cdot \left(\dot{\overrightarrow{\kappa }}-\overrightarrow{\omega }\times \overrightarrow{\kappa }\right)\right]\left(\dot{\overrightarrow{\kappa }}-\overrightarrow{\omega }\times \overrightarrow{\kappa }\right)\right].\end{array}$                            (3.7)

где $T_{VV}$ – коэффициент, характеризующий степень кубической нелинейности диссипативных свойств материла вала.

Отображение инвариантного уравнения с кубической нелинейностью скоростей деформации в модели Кельвина–Фойхта в матричном виде в неподвижном базисе с учетом (3.5) и (3.6) принимает вид:

 $M\left(\kappa ,\dot{\kappa }\right)=E{I}_R\kappa +E{I}_R{T}_V\left[1+B{\left(\dot{\kappa }-\omega R\kappa \right)}^T\left(\dot{\kappa }-\omega R\kappa \right)\right]\left(\dot{\kappa }-\omega R\kappa \right)$,              (3.8)

где $B=T_{VV}/T_V$.

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

4. Учет кубической нелинейности в законе упругости. Аналогично кубическому демпфированию в законе можно учесть и кубическую нелинейность упругих сил по кривизне с коэффициентом $A$:

 $\begin{array}{l}M\left(\kappa ,\dot{\kappa }\right)=E{I}_R\left[1+A\left({\kappa }^T\kappa \right)\right]\kappa +\\ +E{I}_R{T}_V\left[1+B{\left(\dot{\kappa }-\omega R\kappa \right)}^T\left(\dot{\kappa }-\omega R\kappa \right)\right]\left(\dot{\kappa }-\omega R\kappa \right).\end{array}$                                          (4.1)

где $A$ – коэффициент, характеризующий степень нелинейностей упругости материала вала.

5. Уравнения движения вращающегося вала.  Для записи уравнения движения введем вектор прогибов вала $u\left(z,t\right)={\left\{\begin{array}{cc}{u}_x& u_y\end{array}\right\}}^T$, который связан с вектором поворота сечений $\vartheta \left(z,t\right)={\left\{\begin{array}{cc}{\vartheta }_x& {\vartheta }_y\end{array}\right\}}^T$ и вектором кривизны $\kappa \left(z,t\right)={\left\{\begin{array}{cc}{\kappa }_x& {\kappa }_y\end{array}\right\}}^T$ относительно осей $Ox$, $Oy$ дифференциальными соотношениями:

$\vartheta =R{u}^{\prime },\kappa ={\vartheta }^{\prime }=R{u}^{″}$.

Тогда, пренебрегая инерцией поворота сечений при изгибе, уравнения движения вала с погонной массой $m_R$ при действии линейного внешнего трения с коэффициентом $d_e$ могут быть записаны в следующем виде [14]:

 $\left.\begin{array}{l}{m}_R\ddot{u}\left(z,t\right)=Q^{\prime }\left(z,t\right)-d_e\dot{u}\left(z,t\right)\\ 0=M^{\prime }\left(z,t\right)+RQ\left(z,t\right)\end{array}\right\}\Rightarrow m_R\ddot{u}\left(z,t\right)=R{M}^{″}\left(z,t\right)-d_e\dot{u}\left(z,t\right)$,             (5.1)

где $Q\left(z,t\right)$ и $M\left(z,t\right)$ – векторы сосредоточенных поперечных сил и изгибающих моментов соответственно.

Объединяя уравнения (4.1) и (5.1), получим:

 $\begin{array}{l}E{I}_R{{{u^{\prime }}^{\prime }}^{\prime }}^{\prime }+T_{V}E{I}_R\left({{{{\dot{u}}^{\prime }}^{\prime }}^{\prime }}^{\prime }-\omega R{{{u^{\prime }}^{\prime }}^{\prime }}^{\prime }\right)=\\ =-m_R\ddot{u}-E{I}_R{\left[A{G}_{NL}{u}^{″}+T_{VV}{F}_{NL}\left(\dot{u^{″}}-\omega R{u}^{″}\right)\right]}^{\prime \text{}\prime }-d_e\dot{u};\\ F_{NL}=\left({\left(\dot{u^{″}}\right)}^T\dot{u^{″}}-2\omega {\left(\dot{u^{″}}\right)}^{T}R{u}^{″}+{\omega }^2{\left(u^{″}\right)}^T{u}^{″}\right);\\ G_{NL}=\left({\left(u^{″}\right)}^T{u}^{″}\right).\end{array}$                              (5.2)

Переход к безразмерным переменным и величинам осуществляется выбором двух масштабов $U_{\ast }:=l$, $T_{\ast }=\sqrt{m_R{l}^4/E{I}_R}$ и безразмерных комплексов:

$\begin{array}{l}\zeta =z/U_{*},\xi =u/U_{*},\tau =t/T_{*},{\eta }_V=T_V/\left(2T_{*}\right),{\eta }_e=d_e{l}^4/\left(2E{I}_R{T}_{*}\right),\Omega =\omega T_V;\\ {\eta }_{VV}=T_{VV}/\left(2T_{*}{}^3{l}^2\right),\alpha =A/l^2,f_{NL}=F_{NL}{T}_{*}{}^2{l}^2,g_{NL}=G_{NL}{l}^2;\\ \dot{x}=\frac{T_{*}\partial x}{\partial t}=\frac{\partial x}{\partial \tau },x^{\prime }=\frac{U_{*}\partial x}{\partial z}=\frac{\partial x}{\partial \zeta }.\end{array}$ (5.3)

С учетом (5.3) уравнение движения принимает следующий безразмерный вид:

$\begin{array}{l}{{\xi }^{″}}^{″}+2{\eta }_V\left({{\dot{\xi }}^{″}}^{″}-\Omega R{{\xi }^{″}}^{″}\right)=-2{\eta }_e\dot{\xi }-\ddot{\xi }-{\left[\alpha g_{NL}{\xi }^{″}+2{\eta }_{VV}{f}_{NL}\left({\dot{\xi }}^{″}-\Omega R{\xi }^{″}\right)\right]}^{\prime \text{}\prime };\\ f_{NL}=\left({\left({\dot{\xi }}^{″}\right)}^T{\dot{\xi }}^{″}-2\Omega {\left({\dot{\xi }}^{″}\right)}^{T}R{\xi }^{″}+{\Omega }^2{\left({\xi }^{″}\right)}^T{\xi }^{″}\right);\\ g_{NL}={\left({\xi }^{″}\right)}^T{\xi }^{″}.\end{array}$ (5.4)

6. Получение разрешающего уравнения движения. Для сведения уравнения в частных производных к обыкновенному дифференциальному уравнению представим решение в интегральном виде с использованием функции Грина для статического прогиба стержня Бернулли–Эйлера $G_{40}\left(\zeta ,s\right)=\delta \left(\zeta -s\right)$, где $G_{kl}={\partial }^k{\partial }^{l}G\left(\zeta ,s\right)/\partial {\zeta }^k\partial s^l$:

    $\begin{array}{l}\xi \left(\zeta ,\tau \right)+2{\eta }_V\left(\dot{\xi }\left(\zeta ,\tau \right)-\Omega R\xi \left(\zeta ,\tau \right)\right)=\\ ={\displaystyle \underset{0}{\overset{1}{\int }}{G}_{00}\left(\zeta ,s\right)\left[-2{\eta }_e\dot{\xi }\left(s,\tau \right)-\ddot{\xi }\left(s,\tau \right)\right]ds}-\\ -{\displaystyle \underset{0}{\overset{1}{\int }}{G}_{00}\left(\zeta ,s\right){\left[\alpha g_{NL}{\xi }^{″}\left(s,\tau \right)+2{\eta }_{VV}{f}_{NL}\left({\dot{\xi }}^{″}\left(s,\tau \right)-\Omega R{\xi }^{″}\left(s,\tau \right)\right)\right]}^{\prime \text{}\prime }ds}.\end{array}$      (6.1)

Статическая функция Грина для стержня, жестко закрепленного в концевых сечениях, имеет вид:

 $\begin{array}{l}{G}_{00}\left(\zeta ,s\right)=\frac{1}{6}{\left(\zeta -s\right)}^3{\rm H}\left(\zeta -s\right)+\frac{1}{6}\left(-1+3s^2-2s^3\right){\zeta }^3+\frac{1}{2}\left(s-2s^2+s^3\right){\zeta }^2;\\ G_{00}\left(\mathrm{0,}s\right)=\mathrm{0,}{G}_{10}\left(\mathrm{0,}s\right)=\mathrm{0,}{G}_{00}\left(\mathrm{1,}s\right)=\mathrm{0,}{G}_{10}\left(\mathrm{1,}s\right)=\mathrm{0,}\end{array}$

где ${\rm H}\left(\zeta -s\right)$ – функция Хевисайда.

Упрощая выражение (6.1), получаем:

 $\begin{array}{l}\xi \left(\zeta ,\tau \right)+2{\eta }_V\left(\dot{\xi }\left(\zeta ,\tau \right)-\Omega R\xi \left(\zeta ,\tau \right)\right)=\\ ={\displaystyle \underset{0}{\overset{1}{\int }}{G}_{00}\left(\zeta ,s\right)\left[-2{\eta }_e\dot{\xi }\left(s,\tau \right)-\ddot{\xi }\left(s,\tau \right)\right]ds}-\\ -{\displaystyle \underset{0}{\overset{1}{\int }}{G}_{02}\left(\zeta ,s\right)\left[\alpha g_{NL}{\xi }^{″}\left(s,\tau \right)+2{\eta }_{VV}{f}_{NL}\left({\dot{\xi }}^{″}\left(s,\tau \right)-\Omega R{\xi }^{″}\left(s,\tau \right)\right)\right]ds};\\ f_{NL}=\left({\left({\dot{\xi }}^{″}\right)}^T{\dot{\xi }}^{″}-2\Omega {\left({\dot{\xi }}^{″}\right)}^{T}R{\xi }^{″}+{\Omega }^2{\left({\xi }^{″}\right)}^T{\xi }^{″}\right);\\ g_{NL}={\left({\xi }^{″}\right)}^T{\xi }^{″}.\end{array}$               (6.2)

Для записи уравнения (6.2) в матричной форме входящие в него интегралы вычислим приближенно по методу трапеций, разбивая вал на $m$ одинаковых элементов с длиной $h=1/m$ и координатами узлов ${\zeta }_i=\left(i-1\right)h;i=\overline{\mathrm{1,}m+1}$. Значения интеграла в каждом $j$ -ом узле принимается равным

$\begin{array}{l}{\displaystyle \underset{0}{\overset{1}{\int }}G\left({\zeta }_j,s\right)x\left(\tau ,s\right)}ds\approx \frac{h}{2}\left[G_{\left(j\right)\left(1\right)}{x}_{\left(1\right)}+G_{\left(j\right)\left(m+1\right)}{x}_{\left(m+1\right)}\right]+h{\displaystyle \sum _{k=2}^{k=m}{G}_{\left(j\right)\left(k\right)}{x}_{\left(k\right)}},\\ G_{\left(j\right)\left(k\right)}\doteq G\left({\zeta }_j,s_k\right),x_{\left(k\right)}\doteq x\left(\tau ,s_k\right),j,k=\overline{\mathrm{1,}m+1}.\end{array}$

Для того чтобы система не была переопределена, краевые узлы исключаются из рассмотрения, поскольку перемещения в них заранее известны – они нулевые. Таким образом, размерность задачи станет равной $2\left(m-1\right)$.

Введем матрицу функций Грина $G^h{}_{kl}$, вектор узловых перемещений $x$ единичные матрицы $E$ и $I$ размерностью $R^2$ и $R^{m-1}$ соответственно:

 $\begin{array}{l}\underset{\left(m-1\right)\times \left(m-1\right)}{\underbracket{G^h{}_{kl}}}=h\left[\begin{array}{cccc}{G}_{kl}\left({\zeta }_2,s_2\right)& G_{kl}\left({\zeta }_2,s_3\right)& \cdots & G_{kl}\left({\zeta }_2,s_m\right)\\ G_{kl}\left({\zeta }_3,s_2\right)& \ddots & \ddots & ⋮\\ ⋮& \ddots & \ddots & ⋮\\ G_{kl}\left({\zeta }_m,s_2\right)& \cdots & \cdots & G_{kl}\left({\zeta }_m,s_m\right)\end{array}\right],\\ \underset{2\left(m-1\right)\times 1}{x}=\left\{\begin{array}{c}\xi \left({\zeta }_2,\tau \right)\\ \xi \left({\zeta }_3,\tau \right)\\ ⋮\\ \xi \left({\zeta }_m,\tau \right)\end{array}\right\}.\end{array}$

Квадратные матрицы $G^h{}_{kl}$, $I$, $D$ имеют размерность, равную $\left(m-1\right)$, т.е. количеству узлов. Поскольку в каждом узле мы имеем две равнозначные степени свободы – перемещения вдоль осей $x$ и $y$ для матричной записи, каж­дая компонента этих матриц должна быть умножена на единичную матрицу $E$. Для этой операции воспользуемся произведением Кронекера, которое обозначим как $A\otimes B$.

В итоге уравнение в матричной форме принимает следующий вид:

 $\begin{array}{l}{A}_2\ddot{x}+A_1\dot{x}+A_{0}x=F,\\ A_2=G^h{}_{00}\otimes E,\\ A_1=2{\eta }_V\left(I\otimes E\right)+2{\eta }_e\left(G^h{}_{00}\otimes E\right),\\ A_0=I\otimes \left(E-2{\eta }_V\Omega R\right),\\ F=2{\eta }_{VV}\left[{\left(G^h{}_{02}\right)}_f\otimes E\right]\left[-\left(D\otimes E\right)\dot{x}+\Omega \left(D\otimes R\right)x\right]-\\ -\alpha \left[{\left(G^h{}_{02}\right)}_g\otimes E\right]\left[\left(D\otimes E\right)x\right].\end{array}$                                   (6.3)

где ${\left(G^h{}_{02}\right)}_f$ и ${\left(G^h{}_{02}\right)}_g$ – матрицы, полученные умножением столбцов матрицы Грина $G^h{}_{02}$ на соответствующие значения $f_{NL}$ и $g_{NL}$.

Полная нелинейная система уравнений (6.3) зависит от шести параметров $\left\{m,{\eta }_e,{\eta }_V,\Omega ,{\eta }_{VV},\alpha \right\}$, при этом линейная часть – от четырех параметров $\left\{m,{\eta }_e,{\eta }_V,\Omega \right\}$.

Для численного решения сведем уравнение (6.3) к форме Коши:

$\begin{array}{l}\dot{Z}=AZ+f;\\ Z=\left\{\begin{array}{c}x\\ \dot{x}\end{array}\right\},A=\left[\begin{array}{cc}0& E\\ -A_2{}^{-1}{A}_0& -A_2{}^{-1}{A}_1\end{array}\right],f=\left\{\begin{array}{c}0\\ A_2{}^{-1}F\end{array}\right\}.\end{array}$

7. Определение критической скорости. Линейная динамическая система описывается следующим уравнением:

 $A_2\ddot{x}+A_1\dot{x}+A_{0}x=0,x\in R^{2\left(m-1\right)}$         (7.1)

Собственные числа ${\lambda }_k,k=\overline{\mathrm{1,}\mathrm{2,}\dots ,4\left(m-1\right)}$ системы (7.1) являются корнями характеристического уравнения

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

9. Поведение нелинейной динамической системы. Рассмотрим влияние каж­дой из нелинейностей $\left({\eta }_{VV},\alpha \right)$ модели внутреннего демпфирования Кельвина–Фойхта на амплитуды колебаний вала при его вращении с закритической скоростью $\Omega =33>{\Omega }_{\text{crit}}$. В качестве начального возмущения примем перемещения шестого (срединного) узла ${\xi }_x\left({\zeta }_6\mathrm{,0}\right)={\xi }_y\left({\zeta }_6\mathrm{,0}\right)=0.2$. На рис. 3,а показаны перемещения шестого (срединного) узла при значении коэффициента ${\eta }_{VV}={10}^{-5}$ и нулевом коэффициенте нелинейности упругих сил ( $\alpha =0$ ). В этом случае с увеличением числа оборотов вала ${\rm N}$ амплитуды колебаний неограниченно возрастают, стремясь в бесконечность. При учете обеих нелинейностей ( ${\eta }_{VV}={10}^{-5}$, $\alpha =0.015$ ) амплитуда колебаний сначала возрастает и, затем, после, примерно, двадцати оборотов амплитуда устанавливается на постоянное значение, достигнув пятикратного увеличения по сравнению с начальным возмущением (рис. 3,b). Наличие только нелинейности в силах упругости, в данном случае – кубической ( $\alpha =0.015$, ${\eta }_{VV}=0$ ), приводит к ограниченным устойчивым периодическим решениям (рис. 3,c).

Рис. 3. Перемещения срединного узла при закритической скорости при $m=12$, ${\eta }_e=0.02$, ${\eta }_V={10}^{-4}$: (a) ${\eta }_{VV}={10}^{-5}$, $\alpha =0$; (b) ${\eta }_{VV}={10}^{-5}$, $\alpha =0.015$; (c) ${\eta }_{VV}=0$, $\alpha =0.015$.

10. Явление гистерезиса угловой скорости вращения при малом значении времени релаксации неупругих деформаций ${\eta }_V$. Как было показано в предыдущем пункте, наличие кубической нелинейности в силах упругости приводит к появлению дополнительного нетривиального устойчивого решения в закритической области. На рис. 4 приведена бифуркационная диаграмма значений установившихся амплитуд перемещений центрального узла ${\xi }_x\left({\zeta }_6,\tau \right)$ от угловой скорости $\Omega$, которая получена методом установления при следующих значениях параметров: $m=12$, ${\eta }_e=0.02$, ${\eta }_V={10}^{-4}$, ${\eta }_{VV}={10}^{-5}$, $\alpha =0.015$. Видно, что в докритической области, в диапазоне скоростей ${\Omega }_{-}<\Omega <{\Omega }_{\text{crit}}$ существует два устойчивых режима: устойчивое нулевое решение и устойчивое периодическое решение. Эти решения имеют разные притягивающие множества. Можно заключить, что наблюдается субкритическая бифуркация (точка B на рис. 4), неустойчивая ветвь показана схематично красной штриховой линией. Точка А является предельной точкой перехода с верхней устойчивой ветви на нижнюю устойчивую. При угловой скорости выше критической $\Omega >{\Omega }_{\text{crit}}$ существует только одно устойчивое периодическое решение, продолжающее верхнюю ветвь на интервале ${\Omega }_{-}<\Omega <{\Omega }_{\text{crit}}$.

Рис. 4. Бифуркационная диаграмма вдоль верхней периодической ветви, А – предельная точка, B – субкритическая бифуркация ( $m=12$, ${\eta }_e=0.02$, ${\eta }_V={10}^{-4}$, ${\eta }_{VV}={10}^{-5}$, $\alpha =0.015$).

11. Характер прецессии на устойчивой ветви периодических решений. Для определения характера прецессии используется понятие о коэффициенте мгновенной прецессии $\Lambda$, который определяется через скорость изменения полярного угла срединного узла $\alpha \left(\tau \right)$, показанного на рис. 5:

$\dot{\alpha }\left(\tau \right)=\frac{{\dot{\xi }}_y\left({\zeta }_6,\tau \right){\xi }_x\left({\zeta }_6,\tau \right)-{\dot{\xi }}_x\left({\zeta }_6,\tau \right){\xi }_y\left({\zeta }_6,\tau \right)}{{\xi }_x{}^2\left({\zeta }_6,\tau \right)+{\xi }_y{}^2\left({\zeta }_6,\tau \right)}$.

Рис. 5. Прецессия сечения стержня, соответствующего срединному узлу.

Тогда коэффициент мгновенной прецессии $\Lambda$ и его среднее значение $\overline{\Lambda }$ определяются следующим образом:

$\Lambda =\frac{\dot{\alpha }}{\Omega },\overline{\Lambda }=\frac{1}{{\tau }_2-{\tau }_1}{\displaystyle \underset{{\tau }_1}{\overset{{\tau }_2}{\int }}\Lambda \left(s\right)ds}$.

Можно различить следующие случаи: $\overline{\Lambda }>0$ – прямая прецессия, $\overline{\Lambda }<0$ – обратная прецессия, $\overline{\Lambda }=1$ и $\overline{\Lambda }=-1$ – прямая и обратная синхронная прецессии соответственно.

На рис. 6 показано распределение среднего коэффициента прецессии $\overline{\Lambda }$ вдоль верхней устойчивой периодической ветви, в интервале $\Omega >{\Omega }_{-}$ (см. рис. 4) реализуется прямая несинхронная прецессия $0<\overline{\Lambda }<1$, причем скорость вращения вала несколько больше скорости прецессии.

Рис. 6. Распределение среднего значения коэффициента прецессия вдоль верхней устойчивой ветви ( $m=12$, ${\eta }_e=0.02$, ${\eta }_V={10}^{-4}$, ${\eta }_{VV}={10}^{-5}$, $\alpha =0.015$).

На рис. 7 показаны характерные для верхней ветви траектории узлов коллокации с формой изогнутой оси.

Рис. 7. Траектории узлов коллокации при установившемся движении (черные линии) и формы изо-гнутой оси (красные линии) ( $m=12$, ${\eta }_e=0.02$, ${\eta }_V={10}^{-4}$, ${\eta }_{VV}={10}^{-5}$, $\alpha =0.015$, $\Omega =32$).

12. Характер прецессии на устойчивой ветви периодических решений. На рис. 8а показаны значения ${\Omega }_{-}$ и ${\Omega }_{\text{crit}}$ в зависимости от коэффициента внутреннего линейного демпфирования ${\eta }_V$, из которого видно, что при малых значениях ${\eta }_V$ в области ${\Omega }_{-}<\Omega <{\Omega }_{\text{crit}}$ присутствует множественность решений (рис. 8,c), что соответствует субкритической бифуркации. При увеличении коэффициента ${\eta }_V$ значение ${\Omega }_{-}$ постепенно сливается с ${\Omega }_{\text{crit}}$. При некотором значении ${\eta }_V$ складка исчезает: $\Omega ={\Omega }_{\text{crit}}$ – суперкритическая бифуркация (рис. 8,b). Появляется мягкое развитие амплитуд устойчивых периодических (автоколебательных) движений.

Рис. 8. Зависимость ${\Omega }_{-}$ и ${\Omega }_{\text{crit}}$ от коэффициента ${\eta }_V$ (а); распределение амплитуд перемещений срединного узла при субкритической (b) и суперкритической (c) бифуркациях ( $m=12$, ${\eta }_e=0.02$, ${\eta }_{VV}={10}^{-5}$, $\alpha =0.015$).

13. Заключение. В работе показано влияние кубического члена в модели Кельвина–Фойхта для внутреннего демпфирования и кубического члена в законе упругости на динамику вала в докритической и закритической областях скоростей вращения. Без учета нелинейных членов в законе упругости в закритической области прямолинейная форма вала всегда неустойчива (при $\Omega >{\Omega }_{\text{crit}}$: $‖\xi ‖\to \infty$ ). Наличие кубической нелинейности в упругой составляю­щей в деформациях приводит к появлению дополнительной верхней ветви периодических решений в докритической области в диапазоне скоростей $\left[{\Omega }_{-};{\Omega }_{\text{crit}}\right]$, которая продолжается в закритической области. В зависимости от значений времени релаксации неупругих деформаций ${\eta }_V$ в материале вала возможны два различных сценария бифуркаций скорости вращения вала. При относительно большом времени релаксации (здесь ${\eta }_V\ge 0.0004$, рис. 8а) наб­людается сверхкритическая бифуркация прямолинейного вращения вала при критической скорости вращения $\Omega ={\Omega }_{\text{crit}}$. При увеличении скорости вращения происходит мягкое возбуждение асинхронного прецессирования оси вала (рис. 8б) – $0<\overline{\Lambda }<1$. При уменьшении времени релаксации (здесь для ${\eta }_V<0.0004$ рис. 8а) появляется новая точка бифуркации скорости вращения вала ${\Omega }_{-}$ при продолжающемся увеличении критической скорости ${\Omega }_{\text{crit}}$. В диа­пазоне $\left[{\Omega }_{-};{\Omega }_{\text{crit}}\right]$ существует дополнительная ветвь неустойчивых периодических движений (рис. 8,c), соединяющая точки ${\Omega }_{-}$ и ${\Omega }_{\text{crit}}$. Нижняя точка ${\Omega }_{-}$ является бифуркацией типа предельной точкой слияния устойчивого и неустойчивого периодического движения и субкритической бифуркации в точке ${\Omega }_{\text{crit}}$, при достижении которой происходит жесткое возбуждение прецессионного движения. Т.е. возникает область гистерезисного поведения вала: при адиабатическом изменении скорости вращения вала: при движении “вперед–назад” динамическая система проходит через различные состояния.

Исследование выполнено за счет гранта РНФ (проект № 24-19-00333, https://rscf.ru/project/24-19-00333/).