单圆盘转子的临界转速和不平衡响应

转子动力学，钟一鄂
我水平太差，不太好理解这本书，/(ㄒoㄒ)/~~

第一章 http://book.ucdrs.superlib.net/views/specific/2929/bookDetail.jsp?dxNumber=000001034999&d=CE35F309CB85878C55BA563D0038FEED&fenlei=1811050303

第一章 单圆盘转子的临界转速和不平衡响应

刚性支承：轴承与轴承座完全刚性
简单旋转机械可简化为：无质量的弹性转轴上装有圆盘，两端刚性支撑

1.1 转子的涡动

为了避开静变形，将单圆盘转子竖直放置，圆盘位于中间，作转速$\Omega$的匀速运动。此时给转轴加一横向冲击，研究这种运动的性质。
定义坐标系：固定坐标系$xys$，轴向为$s$，其余两个方向为$x$与$y$。圆盘圆心为$o^{\prime }$，静止时圆心所在位置为$o$，$o^{\prime }$的坐标为$(x,y)$。

圆盘此时只受转轴给的弹性恢复力$F=-kr,r=\overline{o{o}^{\prime }}$
圆盘运动的微分方程
$$\{\begin{array}{r}m{x}^{\prime \prime }=-kx\\ m{y}^{\prime \prime }=-ky\end{array}$$
令$w_n^2=k/m$
$$\{\begin{array}{r}{x}^{\prime \prime }+w_n^{2}x=0\\ y^{\prime \prime }+w_n^{2}y=0\end{array}$$
其解可写作
$$\{\begin{array}{r}x=Xcos(w_{n}t+a_x)\\ y=Ysin(w_{n}t+a_y)\end{array}$$
振幅与初相由一瞬间的冲击来决定。
一般情况下，振幅$X\ne Y$，即o’的运动轨迹为椭圆，$o^{\prime }$的运动称为涡动或进动。自然频率为$o^{\prime }$的进动角速度。

复变量$z=x+iy$，用一个复变量代表两个实变量
上面的公式可简化为$z^{\prime \prime }+w_n^{2}z=0$，解为$z=B_1{e}^{i{w}_{n}t}+B_2{e}^{-i{w}_{n}t}$，$B_1,B_2$为复数由一瞬间的冲击来决定。
解的第一项为正进动，运动轨迹为圆，半径为$|B_1|$，运动方向同$\Omega$；解的第一项为反进动，运动轨迹为圆，半径为$|B_2|$，运动方向与$\Omega$相反。两种圆运动合成了椭圆运动。
$o^{\prime }$的运动必然会有下面四种情况：
$B_1=B_2$，轨迹为直线
$B_1\ne B_2$，轨迹为椭圆，正进动半径大于反进动半径时，作正进动，反之，作反进动。（正进动与反进动合成后的运动也可以称为进动）
$B_1=0,B_2\ne 0$，轨迹为圆，反进动
$B_2=0,B_1\ne 0$，轨迹为圆，正进动

$o^{\prime }$的进动属于自然振动，频率就是转轴静止时弯曲振动的自然频率。

如果考虑空气阻力，$z^{\prime \prime }+2n{z}^{\prime }+w_n^{2}z=0$，解为$z=e^{-nt(B_1{e}^{i{w}_n^{\prime }t}+B_2{e}^{-i{w}_n^{\prime }t})},w_n^{\prime }=\sqrt{w_n^2-n^2}$
这种情况下，涡动是衰减的，$o^{\prime }$最终趋于$o$

1.2 圆盘的偏心质量引起的振动，临界转速

圆盘重心$c$与$o^{\prime }$不重合，$o^{\prime }c$为偏心距
$$\{\begin{array}{r}{x}_c^{\prime \prime }=x^{\prime \prime }-e{\Omega }^{2}cos(\Omega t)\\ y_c^{\prime \prime }=y^{\prime \prime }-e{\Omega }^{2}sin(\Omega t)\end{array}$$
由质心运动定理，有
$$\{\begin{array}{r}m{x}_c^{\prime \prime }=-kx\\ m{y}_c^{\prime \prime }=-ky\end{array}$$
可得$o^{\prime }$运动微分方程（强迫振动）
$$\{\begin{array}{r}{x}^{\prime \prime }+w_n^{2}x=e{\Omega }^{2}cos(\Omega t)\\ y^{\prime \prime }+w_n^{2}y=e{\Omega }^{2}sin(\Omega t)\end{array}$$
复变量形式$z^{\prime \prime }+w_n^{2}z=e{\Omega }^2{e}^{i\Omega t}$，特解为$z=A{e}^{i\Omega t}$
将特解带入方程可得振幅$|A|=|\frac{e(\Omega /w_n{)}^2}{1-(\Omega /w_n{)}^2}|$，$z=\frac{e(\Omega /w_n{)}^2}{1-(\Omega /w_n{)}^2}{e}^{i\Omega t}$
可知响应频率与激励频率相同，响应相位与激励相位相同（$\Omega <w_n$）或相差180°（$\Omega >w_n$）
转动过程中$o{o}^{\prime }c$位于同一直线上，直线绕$o$以角速度$\Omega$转动，$o^{\prime }$与$c$作同步正进动
$\Omega <w_n$，$A>0$，$o^{\prime }$与$c$在$o$的同一侧；$\Omega >w_n$，$A<0$，$c$在$o^{\prime }$与$o$的之间，当$\Omega >>w_n$，$A\approx -e$，此时$c$几乎位于$o$，称为自动对心。
$\Omega =w_n$，$A=\infty$，由于公式无阻尼，幅值无限大，$w_n$称为临界角速度，单位为rpm，则称为临界转速。
转子工作转速小于临界转速称为刚性轴，反之称为柔性轴。

如果考虑空气阻力，$z^{\prime \prime }+2n{z}^{\prime }+w_n^{2}z=e{\Omega }^2{e}^{i\Omega t}$，特解为$z=|A|e^{i(\Omega t-\theta )}$
将特解带入微分方程，可求出$|A|$与$\theta$
$|A|=\frac{e(\Omega /w_n{)}^2}{\sqrt{[1-(\Omega /w_n{)}^2{]}^2+(2n/w_n{)}^2(\Omega /w_n{)}^2}}$
$tg(\theta )=\frac{(2n/w_n)(\Omega /w_n)}{1-(\Omega /w_n{)}^2}$
频响曲线为幅值$|A|$与相位$\theta$随频率比$\Omega /w_n$的变化
由于外阻尼，$\Omega /w_n=1$时的响应并不是最大值，最大值发生在$\Omega /w_n\le 1$
在实际中通过测量升速与降速时转子的响应来确定临界转速，因此升速时根据最大响应确定的临界转速比真实的临界转速大，而降速时根据最大响应确定的临界转速比真实的临界转速小。
由于阻尼的存在，$o^{\prime }co$三点并不在同一条直线上。但当$\Omega >>w_n$，$\theta \approx \pi$，此时$c$仍然可认为几乎位于$o$，自动对心。

1.3 圆盘相对于转动坐标系的运动

转子相对转动坐标系的运动为相对运动，转子相对固定坐标系的运动为绝对运动
固定坐标系$sxy$，转动坐标系$s\xi \eta$，以转子角速度$\Omega$绕$s$轴转动
$o^{\prime }$的绝对坐标$z=x+iy$，相对坐标$\zeta =\xi +i\eta$
$z=r{e}^{i(\theta +\Omega t)}$
$\zeta =r{e}^{i\theta }$
因此可以得出 $z=\zeta e^{i\Omega t}$，相对坐标与绝对坐标的关系
带入之前的冲击振动微分方程可得（绝对坐标用相对坐标表示），${\zeta }^{\prime \prime }+2i\Omega {\zeta }^{\prime }+(w_n^2-{\Omega }^2)\zeta =0$，解为$\zeta =B_1{e}^{i(w_n-\Omega )t}+B_2{e}^{-i(w_n+\Omega )t}$
转动坐标系下的解仍是正反进动的合成，但两种进动角速度不同，合成的轨迹也不是椭圆，而是花瓣形。

带入之前的不平衡强迫振动微分方程可得，${\zeta }^{\prime \prime }+2i\Omega {\zeta }^{\prime }+(w_n^2-{\Omega }^2)\zeta =e{\Omega }^2$，解为$\zeta =\frac{e(\Omega /w_n{)}^2}{1-(\Omega /w_n{)}^2}$
$\zeta$与时间无关，即$o^{\prime }$相对转动坐标系不动，而且$|\zeta |=|A|$，固定坐标系下的响应振幅等于转动坐标系下相对平衡位置的大小

1.4 陀螺力矩

圆盘不装在转轴中间，当转轴弯曲，圆盘轴线会与两支点连线有一夹角$\psi$
圆盘对质心$o^{\prime }$的动量矩为$H=J_p\Omega$，动力矩与两支点连线夹角$\psi$
由于进动，圆盘的动量矩方向会不断变化，根据动量矩定理，动量矩变化是受到了力矩（动量变化是受到了力）。
力矩$M_g=-(w_{n}xH)=Hx{w}_n$，力矩方向与平面$o^{\prime }AB$垂直，称为陀螺力矩或回转力矩，是圆盘加于转轴的力矩。

注：一些书上说固定坐标系下叫陀螺力矩，转动坐标系下叫回转力矩。这本书上管这个力矩叫惯性力矩，但我认为在转动坐标系下才叫惯性力矩，地面上的固定坐标系一般认为是惯性坐标系怎么会有惯性力呢？

$|M_g|=H{w}_{n}sin(\psi )$，由于$\psi$很小，$|M_g|=H{w}_n\psi$，力矩与转角成正比，相当于弹性力矩
在正进动$0<\psi <\pi /2$时，陀螺力矩使转轴变形减小，提高了临界角速度
在反进动$\pi >\psi >\pi /2$时，陀螺力矩使转轴变形减大，降低了临界角速度

1.4.1 圆盘的角速度

圆盘不装在转轴中间时，圆盘有绕圆盘直径的转动，也有绕圆盘轴线的转动，用欧拉角来描述。
移动坐标系$o^{\prime }xyz$
固结于圆盘的坐标系$o^{\prime }\zeta \eta \xi$，$o^{\prime }\xi$为圆盘中心轴

圆盘初始位于$o^{\prime }{\xi }_0{\eta }_0{\zeta }_0=o^{\prime }xyz$
绕$o^{\prime }y$转${\theta }_y$到达$o^{\prime }{\xi }_1{\eta }_0{\zeta }_1$
绕$o^{\prime }{\xi }_1$转${\theta }_{\xi }$到达$o^{\prime }{\xi }_1{\eta }_1\zeta$
绕$o^{\prime }\zeta$转$\varphi$到达$o^{\prime }\xi \eta \zeta$
圆盘绝对角速度$w={\theta }_{\xi }^{\prime }+{\theta }_y^{\prime }+{\varphi }^{\prime }$
圆盘相对随动坐标系$o^{\prime }{\xi }_1{\eta }_1\zeta$的转动，$w_1={\theta }_{\xi }^{\prime }+{\theta }_y^{\prime }$

$w$沿随动坐标系$o^{\prime }{\xi }_1{\eta }_1\zeta$各轴的分解
$$\{\begin{array}{r}{w}_{{\xi }_1}={\theta }_{\xi }^{\prime }\\ w_{{\eta }_1}={\theta }_y^{\prime }cos({\theta }_{\xi })\\ w_{\zeta }={\varphi }^{\prime }-{\theta }_y^{\prime }sin({\theta }_{\xi })=\Omega -{\theta }_y^{\prime }sin({\theta }_{\xi })\end{array}$$

$w$沿随动坐标系$o^{\prime }\xi \eta \zeta$各轴的分解
$$\{\begin{array}{r}{w}_{\xi }={\theta }_{\xi }^{\prime }cos(\varphi )+{\theta }_y^{\prime }cos({\theta }_{\xi })sin(\varphi )\\ w_{\eta }=-{\theta }_{\xi }^{\prime }sin(\varphi )+{\theta }_y^{\prime }cos({\theta }_{\xi })cos(\varphi )\\ w_{\zeta }=\Omega -{\theta }_y^{\prime }sin({\theta }_{\xi })\end{array}$$

$w$沿随动坐标系$o^{\prime }\zeta$轴与$o^{\prime }{\xi }_1{\eta }_1$平面的分解
$$\{\begin{array}{r}{w}_e=w_{{\xi }_1}{i}_1+w_{{\eta }_1}{j}_1=w_{\xi }i+w_{\eta }j\\ w_{\zeta }=\Omega -{\theta }_y^{\prime }sin({\theta }_{\xi })\end{array}$$

$w1$沿随动坐标系$o^{\prime }{\xi }_1{\eta }_1\zeta$各轴的分解
$$\{\begin{array}{r}{w}_{1{\xi }_1}={\theta }_{\xi }^{\prime }\\ w_{1{\eta }_1}={\theta }_y^{\prime }cos({\theta }_{\xi })\\ w_{1\zeta }=-{\theta }_y^{\prime }sin({\theta }_{\xi })\end{array}$$

三维旋转：欧拉角、四元数、旋转矩阵、轴角之间的转换 https://zhuanlan.zhihu.com/p/45404840
欧拉角可视化工具 http://danceswithcode.net/engineeringnotes/rotations_in_3d/demo3D/rotations_in_3d_tool.html

1.4.2 圆盘的动量矩

只绕一根轴转动时，动量矩$H=Jw$。圆盘的转动很复杂，因此动量矩也很复杂。
书上之前转动坐标系为$s\xi \eta$，现在为$\zeta \xi \eta$
将动量矩在转动坐标系分解，$G=J_{\xi }{w}_{\xi }i+J_{\eta }{w}_{\eta }j+J_{\zeta }{w}_{\zeta }k$，前两个转动惯量为直径或赤道转动惯量$J_d$，最后一个为极转动惯量$J_p$

轴对称圆盘对其中心$o^{\prime }$的动量矩，$G=J_d{w}_e+J_p{w}_{\zeta }k$
动量矩在自转轴上的投影，$G_{\zeta }=J_p(\Omega -{\theta }_y^{\prime }sin({\theta }_{\xi }))k$
$G=J_d{\theta }_{\xi }^{\prime }{i}_1+J_d{\theta }_y^{\prime }cos({\theta }_{\xi })j_1+J_p(\Omega -{\theta }_y^{\prime }sin({\theta }_{\xi }))k$
转轴截面转角为小量，简化为$G=J_p{w}_{\zeta }k=Hk$

1.4.3 圆盘的动能

$T=Gw/2=(G_{{\xi }_1}{w}_{{\xi }_1}+G_{{\eta }_1}{w}_{{\eta }_1}+G_{\zeta }{w}_{\zeta })/2$
$T=[Jd({\theta }_{\xi }^{\prime 2}+{\theta }_y^{\prime 2}cos({\theta }_{\xi }{)}^2)+J_p{\theta }_y^{\prime 2}sin({\theta }_{\xi }{)}^2+J_p{\Omega }^2-2J_p\Omega {\theta }_y^{\prime }sin({\theta }_{\xi })]/2$
转轴截面转角${\theta }_{\xi }$与${\theta }_y$为小量，则$sin({\theta }_{\xi })\approx {\theta }_{\xi }\approx {\theta }_x,cos({\theta }_{\xi })\approx 1$，忽略二阶小量，简化为$T=[Jd({\theta }_x^{\prime 2}+{\theta }_y^{\prime 2})+J_p{\Omega }^2-2J_p\Omega {\theta }_y^{\prime }{\theta }_x]/2$

注：无穷小量Infinitesimals即以数0为极限的变量，无限接近于0，二阶小量是其中一种。

1.5 圆盘绕其中心的转动方程

动量矩定理：动量矩的变化量等于外力矩
移动坐标系$o^{\prime }xyz$中动量矩的变化率$\frac{dG}{dt}=M$
转动坐标系$o^{\prime }{\xi }_1{\eta }_1\zeta$中动量矩的变化率$\frac{d\hat{G}}{dt}=\frac{dG}{dt}-w_1\times G$，$w_1$为转动坐标系的角速度。
$-(w_1\times G)$为陀螺力矩

将转动坐标系中的动量矩变化率沿转动坐标系各轴分解
$$\{\begin{array}{r}\frac{d{G}_{{\xi }_1}}{dt}+w_{1{\eta }_1}{G}_{\zeta }-w_{1\zeta }{G}_{{\eta }_1}=M_{{\xi }_1}\\ \frac{d{G}_{{\eta }_1}}{dt}+w_{1\zeta }{G}_{\xi }-w_{1{\xi }_1}{G}_{\zeta }=M_{{\eta }_1}\\ \frac{d{G}_{\zeta }}{dt}+w_{1{\xi }_1}{G}_{{\eta }_1}-w_{1{\eta }_1}{G}_{{\xi }_1}=M_{\zeta }\end{array}$$
$$\{\begin{array}{r}{J}_d({\theta }_{\xi }^{\prime \prime }+{\theta }_y^{\prime 2}sin({\theta }_{\xi })cos({\theta }_{\xi }))+H{\theta }_y^{\prime }cos({\theta }_{\xi })=M_{{\xi }_1}\\ J_d({\theta }_y^{\prime \prime }cos({\theta }_{\xi })-2{\theta }_{\xi }{\theta }_y^{\prime }sin({\theta }_{\xi }))-H{\theta }_{\xi }^{\prime }=M_{{\eta }_1}\\ \frac{dH}{dt}=M_{\zeta }\end{array}$$
转轴截面转角${\theta }_{\xi }$与${\theta }_y$为小量，则$sin({\theta }_{\xi })\approx {\theta }_{\xi }\approx {\theta }_x,cos({\theta }_{\xi })\approx 1$，忽略二阶小量，简化为
$$\{\begin{array}{r}{J}_d{\theta }_x^{\prime \prime }+H{\theta }_y^{\prime }=M_{{\xi }_1}\approx M_x\\ J_d{\theta }_y^{\prime \prime }-H{\theta }_x^{\prime }=M_{{\eta }_1}\approx M_y\\ \frac{dH}{dt}=M_{\zeta }\end{array}$$
将圆盘的转动微分方程从圆盘的转动坐标系$o^{\prime }{\xi }_1{\eta }_1\zeta$转移到转动坐标系$o^{\prime }xy\zeta$

转子正常运转时，驱动力矩与阻力矩平衡，$M_{\zeta }=0$
$$\{\begin{array}{r}{J}_d{\theta }_x^{\prime \prime }+J_p\Omega {\theta }_y^{\prime }=M_x\\ J_d{\theta }_y^{\prime \prime }-J_p\Omega {\theta }_x^{\prime }=M_y\end{array}$$

圆盘所受的外力与外力矩需要转轴给的弹性力与弹性力矩来平衡，在转动坐标系中，外力为惯性力与惯性力矩
圆盘的运动微分方程，包括圆盘的移动与转动微分方程，圆盘$o^{\prime }$有四个位移，绕x与绕y的转角，沿x与y的位移，
$$\{\begin{array}{r}m{x}^{\prime \prime }+k_{11}x+k_{14}{\theta }_y=0\\ m{y}^{\prime \prime }+k_{22}y-k_{23}{\theta }_x=0\\ J_d{\theta }_x^{\prime \prime }+J_p\Omega {\theta }_y^{\prime }+k_{32}x+k_{33}{\theta }_x=0\\ J_d{\theta }_y^{\prime \prime }-J_p\Omega {\theta }_x^{\prime }+k_{41}x+k_{41}{\theta }_y=0\end{array}(1.47)$$
用拉格朗日方程也可以推出1.47

1.6 考虑陀螺力矩时，转子的临界角速度

求解1.47四个微分方程的特征根可得到转子振动的自然频率$w_n$，即进动角速度。
临界角速度就是与进动角速度相等的(工作)转动角速度。

如果转轴截面为圆，则
$$\{\begin{array}{r}{k}_{11}=k_{22}=k_{rr}\\ k_{33}=k_{44}=k_{\psi \psi }\\ k_{14}=k_{41}=k_{23}=k_{32}=k_{\psi r}=k_{r\psi }\end{array}$$
复变量$z=x+iy,\psi ={\theta }_y-i{\theta }_x$
则1.47变为
$$\{\begin{array}{r}m{z}^{\prime \prime }+k_{rr}z+k_{r\psi }\psi =0\\ J_d{\psi }^{\prime \prime }-i{J}_p\Omega \psi +k_{\psi r}z+k_{\psi \psi }\psi =0\end{array}$$
$$\{\begin{array}{r}{z}^{\prime \prime }+w_{rr}^{2}z+w_{r\psi }^2\psi =0\\ {\psi }^{\prime \prime }-i\frac{J_p}{J_d}\Omega \psi +w_{\psi r}^{2}z+w_{\psi \psi }^2\psi =0\end{array}$$
解为$z=z_0{e}^{i{w}_{n}t},\psi ={\psi }_0{e}^{i{w}_{n}t}$
带入方程可得
$$\{\begin{array}{r}(-w_n^2+w_{rr}^2)z_0+w_{r\psi }^2{\psi }_0=0\\ w_{\psi r}^2{z}_0+[-w_n^2+\frac{J_p}{J_d}\Omega w_n+w_{\psi \psi }^2]{\psi }_0=0\end{array}$$
将上面复变量微分方程组第一个带入第二个中，得到特征值方程（频率方程），微分方程组变成了一元四次方程
$(-w_n^2+w_{rr}^2)[-w_n^2+\frac{J_p}{J_d}\Omega w_n+w_{\psi \psi }^2]-w_{r\psi }^2{w}_{\psi r}^2=0$
该方程有四个根，四个进动角速度

如果没有陀螺力矩
$$\{\begin{array}{r}(-w_n^2+w_{rr}^2)z_0+w_{r\psi }^2{\psi }_0=0\\ (-w_n^2+w_{\psi \psi }^2){\psi }_0+w_{\psi r}^2{z}_0=0\end{array}$$
$(-w_n^2+w_{rr}^2)(-w_n^2+w_{\psi \psi }^2)-w_{r\psi }^2{w}_{\psi r}^2=0$
$w_n^4-(w_{rr}^2+w_{\psi \psi }^2)w_n^2+w_{rr}^2{w}_{\psi \psi }^2-w_{r\psi }^2{w}_{\psi r}^2=0$
该方程有四个根，四个进动角速度，但进动角速度不会变化

进动角速度的振型，$\frac{z}{\psi }=\frac{z_0}{{\psi }_0}=\frac{-w_{r\psi }^2}{w_{rr}^2-w_n^2}$

example

圆盘质量$m=20kg$，半径$R=12cm$；转轴跨度$l=75cm$，直径$d=3cm$。圆盘至左支点的距离$a=l/3=25cm$。求转子的进动角速度，振型，临界角速度。

求解进动角速度与临界角速度

# -*- coding: utf-8 -*-
"""
@time: 2021-08-22 下午 04:51
@author: leslie lee

转轴a位置受P力时 a处的挠度与截面转角
r = 4*P*l**3/(243*E*I)
psi = 2*P*l**2/(81*E*I)
转轴a位置受M力矩时 a处的挠度与截面转角
r = 2*M*l**2/(81*E*I)
psi = M*l/(9*E*I)
如果 P=1 M=1 得到的位移就是柔度系数
刚度矩阵为柔度矩阵的逆

sympy 复数结果 提取实部 ? 没有找到简单的好方法
as_real_imag() 结果为复数
re() 提取实部
im() 提取虚部
np.real()无法处理 应该是无法识别sympy中的虚数单位 I
"""
import numpy as np
import sympy as sp


m = 20 # 圆盘质量 kg
R = 12 # 圆盘半径 cm
l = 75 # 转轴长度 cm
d = 3 # 转轴直径 cm
a = l/3 # 圆盘位置 cm
E = 20.58*1e6  # 转轴弹性模量 N/cm^2  钢

# 圆盘转动惯量
Jp = m*R**2/2
Jd = Jp/2
# 转轴截面惯性矩
I = np.pi*d**4/64
# 转轴
P = 1
a_rr = 4*P*l**3/(243*E*I) # 单位力 产生的位移
a_rp = 2*P*l**2/(81*E*I) # 单位力 产生的转角
M = 1
a_pr = 2*M*l**2/(81*E*I) # 单位力矩 产生的位移
a_pp = M*l/(9*E*I) # 单位力矩 产生的转角
A = np.array([[a_rr, a_rp], [a_pr, a_pp]]) # 柔度矩阵 单位力产生的位移 cm/N rad/N cm/M rad/M
K = np.linalg.pinv(A) # 刚度矩阵 产生单位位移所需的力 N/cm N/rad M/rad M/cm
k_rr, k_rp, k_pr, k_pp = K[0,0], K[0,1], K[1,0], K[1,1]

""" 单位换算
krr N/cm
krp=kpr N  N/rad=N  M/cm=N*cm/cm=N   rad无量纲
kpp N*cm

w2_rr = k_rr/m   (N/cm)/(kg)=(kg*m/s^2/cm)/(kg)=100*1/s^2
w2_rp = k_rp/m   (N)/(kg)=(kg*m/s^2)/(kg)=100*cm/s^2
w2_pp = k_pp/Jd  (N*cm)/(kg*cm^2)=(kg*m/s^2*cm)/(kg*cm^2)=100*1/s^2
w2_pr = k_pr/Jd  (N)/(kg*cm^2)=(kg*m/s^2)/(kg*cm^2)=100*1/(s^2*cm)
"""
w2_rr = k_rr/m*1e2
w2_rp = k_rp/m*1e2
w2_pp = k_pp/Jd*1e2
w2_pr = k_pr/Jd*1e2


# Omega = 200 # 转速 rad/s
wn = sp.symbols('wn') # 特征根 1/s    wn = sp.Symbol('wn')
# 遍历
for Omega in np.arange(0,2000,200):
    # 方程1 有陀螺力矩
    eq1 = (-wn**2 + w2_rr)*(-wn**2 + (Jp/Jd)*Omega*wn + w2_pp) - w2_rp*w2_pr
    # 方程2 无陀螺力矩
    eq2 = (-wn**2 + w2_rr)*(-wn**2 + w2_pp) - w2_rp*w2_pr
    # 解方程

    print(sp.solve(eq2))

求临界转速，我只列举两个方法：
方1：可以绘制坎贝尔图观察
方2：$Omega$为反向，将$w_n=-\Omega$带入方程求解，$Omega$为正向，将$w_n=\Omega$带入方程求解

Omega = sp.symbols('Omega') # 特征根 1/s    Omega = sp.Symbol('Omega')
eq1 = (-Omega**2 + w2_rr)*(-Omega**2 + (Jp/Jd)*Omega*Omega + w2_pp) - w2_rp*w2_pr
eq2 = (-Omega**2 + w2_rr)*(-Omega**2 - (Jp/Jd)*Omega*Omega + w2_pp) - w2_rp*w2_pr

print(sp.solve(eq1))
print(sp.solve(eq2))
'''
[-244.442111078112, 244.442111078112, -1420.23599311756*I, 1420.23599311756*I]
[-844.987991611588, -237.205839436574, 237.205839436574, 844.987991611588]
'''

前三阶临界转速为：反向237.205 正向244.442 反向844.988 rad/s，且只有三阶（取复数的实部）

求解振型

# 角速度为200rad/s时的四阶进动角速度
wns = [243.8324214, 1652.071068, -237.8166687, -1258.086821]
for wn in wns:
    z_psi = -w2_rp/(w2_rr - wn**2)  # 单位 cm
    print(z_psi)

'''
50.93070501720033
-0.5577122225702382
46.29261349959544
-0.9856062550798872
'''

振型图的绘制，以角速度为$200rad/s$为例：
求出$o^{\prime }$处的振型，两端为简支，已知一阶进动与二阶进动的振型形状，则可以绘制出大概的振型图

如果圆盘装在中间，那么$k_{rr}=48EI/l^3,k_{\psi \psi }=12EI/l,k_{r\psi }=k_{\psi r}=0$
微分方程为
$$\{\begin{array}{r}m{z}^{\prime \prime }+k_{rr}z=0\\ J_d{\psi }^{\prime \prime }-iH\psi +k_{\psi \psi }\psi =0\end{array}$$
自然频率为
$$\{\begin{array}{r}{w}_{nz}=\sqrt{k_{rr}/m}\\ -w_{n\psi }^2+\frac{J_p}{J_d}\Omega w_{n\psi }+w_{\psi \psi }^2=0\end{array}$$

import numpy as np
import sympy as sp


m = 20 # 圆盘质量 kg
R = 12 # 圆盘半径 cm
l = 75 # 转轴长度 cm
d = 3 # 转轴直径 cm
a = l/2 # 圆盘位置 cm
E = 20.58*1e6  # 转轴弹性模量 N/cm^2  钢

# 圆盘转动惯量
Jp = m*R**2/2
Jd = Jp/2
# 转轴截面惯性矩
I = np.pi*d**4/64
# 转轴
k_rr, k_rp, k_pr, k_pp = 48*E*I/l**3, 0, 0, 12*E*I/l

w2_rr = k_rr/m*1e2
w2_rp = k_rp/m*1e2
w2_pp = k_pp/Jd*1e2
w2_pr = k_pr/Jd*1e2


Omega = sp.symbols('Omega') # 特征根 1/s    Omega = sp.Symbol('Omega')
eq1 = -Omega**2 + Jp/Jd*Omega*Omega + w2_pp
eq2 = -Omega**2 - Jp/Jd*Omega*Omega + w2_pp

print((w2_rr)**0.5)
print(sp.solve(eq1))
print(sp.solve(eq2))

'''
215.75649075713136
[-1348.47806723207*I, 1348.47806723207*I]
[-778.544175112742, 778.544175112742]
'''

因此如果圆盘置中，可以求出三个自然频率，第一个方程一个，第二个方程两个。第二个微分方程求出的也是两个进动角速度。
临界转速有两个，一阶临界转速为215.756，二阶临界转速为778.544（反进动）

1.7 弹性支承对转子临界转速的影响

轴承考虑弹性，左右支承点也有位移了就不再是之前的0，使得转子的进动角速度或临界转速降低。
支点坐标$A^{\prime }=(x_A,y_A),B^{\prime }=(x_B,y_B)$，$o^{\prime }$的坐标$(x,y)$
$o^{\prime }$的坐标以及截面转角$(x,y,{\theta }_x,{\theta }_y)=(x_1+x^{\prime },y_1+y^{\prime },{\theta }_{xA}+{\theta }_x^{\prime },{\theta }_{yA}+{\theta }_y^{\prime })$
$o^{\prime }$的位移是由两部分叠加的，第一部分为轴承变形造成的转轴位移（刚体位移），第二部分为转轴弯曲变形。

下面考虑当转轴受力而只有支点发生形变时，即受外力后，弹性力来自轴承而非轴承加转轴。求转轴的柔度系数。
支点位移即$o^{\prime }$位移
$$\{\begin{array}{r}-{\theta }_{xA}=(y_B-y_A)/l\\ {\theta }_{yA}=(x_B-x_A)/l\\ x_1=(1-a/l)x_A+(a/l)x_B\\ y_1=(1-a/l)y_A+(a/l)y_B\end{array}$$
支点弹簧沿x与y方向刚度相同时，可以用复变量表示
$$\{\begin{array}{r}{r}_1=x_1+i{y}_1=(1-a/l)r_A+(a/l)r_B\\ {\psi }_A={\theta }_{yA}-i{\theta }_{xA}=(r_B-r_A)/l\end{array}$$
$r_A$与$r_B$为转轴两端简支时求出的支反力与轴承刚度之比，$r_A=R_A/k_A$与$r_B=R_B/k_B$，$k_A$与$k_A$为轴承刚度
受外力时的支反力
$$\{\begin{array}{r}{R}_A=P(l-a)/l\\ R_B=Pa/l\end{array}1.7.1$$
受外力矩时的支反力
$$\{\begin{array}{r}{R}_A=-M/l\\ R_B=M/l\end{array}1.7.2$$
将1.7.1带入$o^{\prime }$位移，得$r_1,{\psi }_A$，$r_1=a_{rr}P,{\psi }_A=a_{\psi r}P$
将1.7.2带入$o^{\prime }$位移，得$r_1,{\psi }_A$，$r_1=a_{r\psi }M,{\psi }_A=a_{\psi \psi }M$
这是只有支点位移时的柔度系数（只有轴承变形），再加上支点0位移时的柔度系数（只有转轴变形），即为总柔度系数。

以之前的实例，研究轴承刚度变化对临界转速的影响
将之前求得的转轴柔度系数与支点位移时的柔度系数叠加，得到总的柔度系数
由柔度矩阵取逆得出刚度矩阵，将刚度系数带入特征方程可求出自然频率与临界角速度

import numpy as np
import sympy as sp


m = 20 # 圆盘质量 kg
R = 12 # 圆盘半径 cm
l = 75 # 转轴长度 cm
d = 3 # 转轴直径 cm
a = l/3 # 圆盘位置 cm
E = 20.58*1e6  # 转轴弹性模量 N/cm^2  钢

# 圆盘转动惯量
Jp = m*R**2/2
Jd = Jp/2
# 转轴截面惯性矩
I = np.pi*d**4/64
# 转轴
kc = 81*E*I/l**3
kA = kB = np.inf # kc/10

# 支点变形时转轴的柔度系数
a_rr1 = 1/kA*(1-a/l)**2 + 1/kB*(a/l)**2
a_pr1 = (1/kB*a/l - 1/kA*(1-a/l))/l
a_rp1 = (1/kB*a/l - 1/kA*(1-a/l))/l
a_pp1 = (1/kA + 1/kB)/l**2
# 转轴变形的柔度系数
a_rr2 = 4*l**3/(243*E*I)
a_rp2 = 2*l**2/(81*E*I)
a_pr2 = 2*l**2/(81*E*I)
a_pp2 = l/(9*E*I)
# 合并
a_rr = a_rr1 + a_rr2
a_pr = a_pr1 + a_pr2
a_rp = a_rp1 + a_rp2
a_pp = a_pp1 + a_pp2

A = np.array([[a_rr, a_rp], [a_pr, a_pp]]) # 柔度矩阵 单位力产生的位移 cm/N rad/N cm/M rad/M
K = np.linalg.pinv(A) # 刚度矩阵 产生单位位移所需的力 N/cm N/rad M/rad M/cm
k_rr, k_rp, k_pr, k_pp = K[0,0], K[0,1], K[1,0], K[1,1]

w2_rr = k_rr/m*1e2
w2_rp = k_rp/m*1e2
w2_pp = k_pp/Jd*1e2
w2_pr = k_pr/Jd*1e2

Omega = sp.symbols('Omega') # 特征根 1/s    Omega = sp.Symbol('Omega')
eq1 = (-Omega**2 + w2_rr)*(-Omega**2 + (Jp/Jd)*Omega*Omega + w2_pp) - w2_rp*w2_pr
eq2 = (-Omega**2 + w2_rr)*(-Omega**2 - (Jp/Jd)*Omega*Omega + w2_pp) - w2_rp*w2_pr

print(sp.solve(eq1))
print(sp.solve(eq2))

可以发现，减小支承刚度会显著的降低临界转速。