转子动力学与模态分析：临界转速计算与稳定状态下的轴心轨迹解析，圆盘质量不平衡条件下的振动响应及...
转子动力学，模态分析求临界转速，稳定状态下转子轴心轨迹，圆盘质量不平衡条件下振动响应+操作视频。
转子系统跑着跑着突然剧烈振动？八成是遇到临界转速了。这玩意儿就像荡秋千到最高点突然有人推你一把，系统固有频率和旋转频率共振时，设备能给你表演原地升天。今天咱们直接开撸代码，手把手带你看转子怎么跳舞。
先整点硬核的，临界转速计算直接上特征值分析。敲个Python脚本感受下：
import numpy as np from scipy.linalg import eig m = 50 # 质量 kg k = 1e6 # 刚度 N/m c = 500 # 阻尼 N·s/m # 构建质量矩阵和刚度矩阵 M = np.array([[m, 0], [0, m]]) K = np.array([[k, -k], [-k, k]]) C = np.array([[c, -c], [-c, c]]) # 特征值求解 omega = 2000 # 初始猜测转速 rad/s A = np.vstack([np.hstack([np.zeros_like(M), np.eye(2)]), np.hstack([-np.linalg.inv(M)@K, -np.linalg.inv(M)@C])]) w, v = eig(A) # 特征值分解 critical_speeds = np.abs(np.imag(w))/(2*np.pi) print(f'临界转速：{critical_speeds[np.argmax(critical_speeds)]:.2f} Hz')这段代码暗藏玄机：把二阶微分方程转化为状态空间方程，特征值的虚部直接对应系统固有频率。注意这里没考虑陀螺效应，实际做转子动力学得把回转力矩矩阵加进去，不然计算结果会漏掉重要模态。
画轴心轨迹才是真的带劲，上实测数据处理代码：
import matplotlib.pyplot as plt from scipy.signal import hilbert # 假设x,y是振动传感器数据 t = np.linspace(0, 10, 1000) x = 0.1*np.sin(2*np.pi*25*t) + 0.05*np.random.randn(len(t)) y = 0.08*np.cos(2*np.pi*25*t) + 0.03*np.random.randn(len(t)) # 希尔伯特变换提取包络 analytic_x = hilbert(x) analytic_y = hilbert(y) plt.figure(figsize=(8,6)) plt.plot(x, y, alpha=0.3, label='原始轨迹') plt.plot(np.real(analytic_x), np.real(analytic_y), 'r--', lw=2, label='包络线') plt.xlabel('X方向位移 (mm)') plt.ylabel('Y方向位移 (mm)') plt.title('轴心轨迹-稳定状态') plt.axis('equal') plt.legend() plt.show()看到没？原始数据乱得像醉酒走位，包络线一画立马现出原形——典型的椭圆轨迹说明存在固定方向的激振力。如果轨迹变成内八字或者外八字，那可得赶紧查是不是碰摩故障了。
质量不平衡响应模拟最刺激，直接上微分方程暴力求解：
from scipy.integrate import solve_ivp def rotor_eq(t, z, m, c, k, e, omega): x, x_dot, y, y_dot = z Fx = m*e*omega**2 * np.cos(omega*t) Fy = m*e*omega**2 * np.sin(omega*t) dxdt = x_dot dydt = y_dot dx_dot = (Fx - c*x_dot - k*x)/m dy_dot = (Fy - c*y_dot - k*y)/m return [dxdt, dx_dot, dydt, dy_dot] # 模拟参数 e = 0.001 # 偏心距1mm res = solve_ivp(rotor_eq, [0, 20], [0,0,0,0], args=(m, c, k, e, 2*np.pi*25), method='LSODA', rtol=1e-5) # 绘制时域波形 plt.plot(res.t, res.y[0], label='X方向') plt.plot(res.t, res.y[2], label='Y方向') plt.title('质量不平衡振动响应') plt.xlabel('时间 (s)') plt.ylabel('振幅 (mm)') plt.legend()这里用龙格库塔法硬刚非线性方程，注意看振幅随时间的变化趋势——刚开始乱颤那是瞬态响应，后面走成稳定正弦波才是稳态。代码里偏心距参数e调大点，分分钟给你画出振幅爆表的刺激曲线。
操作视频要点：在轴心轨迹图上加个实时更新的红点，用matplotlib.animation模块整成动态图。记得设置interval参数控制刷新速度，ffmpeg输出视频时要选对编码格式，不然生成的视频可能只有几KB——别问我是怎么知道的。