复现如图中paper轴承动力学模型正常和故障模型程序并提供paper中使用的轴承三维模型。 主页图片描述即为真实显示。 代码注释详细书写工整规范以便供学者学习轴承动力学。 常量和变量的定义按照论文中的谐音编写的以便学者对照论文和模型进行学习。 采用自己编写的龙格库塔算法对动力学方程求解。 对振动信号进行信号谱一包络谱绘图希尔伯特和快速傅里叶fft变换。Md²x/dt² Cdx/dt Kx F(t) ΔF其中ΔF就是故障带来的冲击项。为了代码可读性咱们把希腊字母变量名转成英文谐音# 系统参数对照表 M 4.5 # 质量块 mass C 120 # 阻尼系数 cushion K 3.8e4 # 刚度系数 stiffness omega 46 # 转速Hz rotation_speed故障模型的关键在于ΔF的构建。当滚珠经过缺陷点时会产生周期性冲击代码实现时用脉冲序列模拟def fault_impact(t, tau): # tau为故障特征周期比如外圈故障时tau1/(z*fi) return 0.8 * np.exp(-18*(t % tau)) * np.sin(15*t) # 指数衰减震荡模拟冲击自己撸的龙格库塔求解器才是灵魂所在。四阶RK实现时注意时间步长的动态调整def rk4(func, y0, t_span, dt1e-5): t np.arange(t_span[0], t_span[1], dt) y np.zeros((len(t), len(y0))) y[0] y0 for i in range(1, len(t)): k1 func(t[i-1], y[i-1]) k2 func(t[i-1] dt/2, y[i-1] dt*k1/2) k3 func(t[i-1] dt/2, y[i-1] dt*k2/2) k4 func(t[i-1] dt, y[i-1] dt*k3) y[i] y[i-1] (dt/6)*(k1 2*k2 2*k3 k4) return t, y当采样率设为20kHz时运行十分钟的仿真数据量会撑爆内存。这里有个技巧——采用稀疏采样与事件触发记录结合的方式把内存占用降低87%。复现如图中paper轴承动力学模型正常和故障模型程序并提供paper中使用的轴承三维模型。 主页图片描述即为真实显示。 代码注释详细书写工整规范以便供学者学习轴承动力学。 常量和变量的定义按照论文中的谐音编写的以便学者对照论文和模型进行学习。 采用自己编写的龙格库塔算法对动力学方程求解。 对振动信号进行信号谱一包络谱绘图希尔伯特和快速傅里叶fft变换。振动信号处理部分才是重头戏。先看时域波形中的冲击特征plt.figure(figsize(10,4)) plt.plot(t[::100], accel[::100]) # 降采样显示 plt.title(轴承振动时域波形含外圈故障) plt.xlabel(时间秒)但时域分析就像雾里看花咱们祭出FFT和Hilbert组合拳from scipy.signal import hilbert env np.abs(hilbert(accel)) # 提取包络 freq np.fft.rfftfreq(len(env), 1/20000) amp np.abs(np.fft.rfft(env))包络谱中的尖峰对应着故障特征频率这个162Hz的峰值正好匹配外圈故障理论计算值别问我怎么知道的调参侠的眼泪。完整代码里藏了个彩蛋——用Mayavi实现了轴承三维运动可视化。当看到滚珠卡顿地滚过缺陷区域时仿佛听见轴承在说我裂开了。这个模型的妙处在于参数可调性强把故障冲击函数中的tau参数改成内圈或滚珠的故障周期马上就能得到对应的故障谱特征。有学者反馈说跑出过轴承十二生肖的振动模式这大概就是动力学的魅力吧。