简介本资源面向机械工程、故障诊断与振动分析领域的高校师生及工业设备维护工程师聚焦轴承动力学建模与早期故障识别这一核心问题提供基于MATLAB的可运行仿真与诊断方案。压缩包共6个文件全部为.m脚本含vdp类与v类命名模型文件总大小仅3KB轻量紧凑涵盖轴承非线性动力学方程构建、ODE45高精度数值求解、振动响应仿真及典型故障特征提取等关键环节代码结构清晰、注释友好便于理解滚动体-套圈接触力学建模逻辑与故障信号演化规律。已有3009人学习下载适用于课程设计、科研建模入门及预测性维护算法验证场景读者可直接调用脚本复现轴承在正常、内圈/外圈缺陷等工况下的时频域动态响应并结合信号处理思路开展峭度分析、谱图判读等诊断实践是贯通理论建模与工程诊断的实用型MATLAB教学与开发参考。 做轴承故障诊断的工程人员和研究生基本都逃不过一个坎手里只有一堆振动数据想验证诊断算法却拿不到带故障的轴承样本。买一套故障轴承不便宜试验台搭建周期又长而且故障程度还不一定可控。我的做法是用MATLAB先做轴承动力学建模用ODE45把运动方程解出来仿真出一组“带故障”的振动信号再拿它调试包络谱、特征频率识别这些诊断逻辑。这套流程我用了很多次稳定、可控而且能重复生成任意工况的数据比纯理论推导直观也比纯数据驱动的方法更容易定位问题。这篇文章就完整记录一下从模型建立、参数设置、ODE45求解到故障诊断代码实现希望给正在做相关方向的朋友提供一个可以直接参考的起点。1. 轴承动力学模型从哪来问题定义与建模思路1.1 轴承故障诊断为什么需要动力学模型直接回答一个问题既然有振动信号为什么不直接做频谱分析非要先建模我的理解是实际采集到的轴承信号里混着转频、齿轮啮合频率、环境噪声甚至传感器安装共振很难干净地把故障特征分离出来。动力学模型的价值在于它把“故障”变成一个可量化的激励源你可以控制故障位置、故障程度、转速和负载然后观察振动响应的变化规律。这个过程中建立的“特征频率-故障类型”对应关系比单纯从数据里翻出来的规律可靠得多也更容易解释。另外在开发诊断算法时真实故障数据是稀缺资源。你不可能为了调一个滤波器参数就去买几十个故障轴承。仿真模型相当于一个“万能试验台”你可以批量生成数据然后测试不同的诊断方法时域指标、包络谱、倒频谱、深度学习分类器都能在这套数据上先跑通。我自己的习惯是先用仿真数据验证算法逻辑再拿实际数据测试这样能省下大量调试时间。1.2 单自由度质量-弹簧-阻尼模型够不够用最早我尝试过直接用有限元模型做一个完整的三维滚动轴承网格一划计算一次要好几个小时而且参数调整极不灵活。后来慢慢发现对于故障诊断这种“看特征频率”的场景单自由度模型已经完全够用。它的物理意义非常清晰把轴承-转子系统简化为一个等效质量m通过一个等效刚度k和阻尼c支承在基础上。当滚动体经过局部缺陷时会产生一个周期性冲击力F(t)系统的振动响应就是对这个冲击的衰减振荡。最基本的运动方程写出来是$$m\ddot{x}c\dot{x}kxF(t)$$其中x可以理解为轴承外圈或传感器安装位置的等效位移。m、c、k不是直接从轴承手册里查到的而是根据系统共振行为反推的等效值。实际操作中先估计系统的固有频率fn和阻尼比zeta然后令km*(2pifn)^2c2zetasqrt(k*m)。这样做的好处是你只需要关注你关心的频带模型动态特性基本可控。对于包络谱分析我们通常关注的是冲击激起的系统高频共振所以固有频率可以设置在1kHz到5kHz之间这也是很多实际诊断中带通滤波器常用的频带。1.3 故障冲击激励怎么模拟才真实故障激励不是纯粹的脉冲函数因为滚动体通过缺陷时接触状态是渐变的。常用的做法是用一个带高频衰减振荡的短时冲击来表示。也就是说每次冲击不是理想冲激而是激发系统固有频率后快速衰减的振荡。如果直接用一个Dirac函数ODE45会因为阶跃过于尖锐而需要极小的步长算得很慢也不符合物理过程。我通常把故障力建模成如下形式$$F_f(t)\sum_{i} A_i \exp\left[-\zeta_f \omega_n (t-t_i)\right]\sin\left(\omega_n \sqrt{1-\zeta_f^2}(t-t_i)\right) \cdot u(t-t_i)$$其中t_i是第i次冲击的发生时刻A_i是冲击幅值zeta_f控制冲击衰减速度omega_n是冲击响应的角频率。这个形式其实模拟了“滚动体撞击缺陷边缘-激发局部共振-逐渐衰减”的过程。冲击时刻t_i按照故障特征频率的倒数间隔产生也就是每经过一个缺陷产生一次冲击。为了让仿真更贴近真实信号我还会在激励上加入小幅度的随机扰动模拟转速波动和载荷变化。2. 用ODE45求解运动方程从数学方程到数值解2.1 二阶方程转一阶状态空间ODE45是MATLAB中基于Dormand-Prince方法的非刚性常微分方程求解器它要求你把待求解的方程写成一阶微分方程组的形式。很多初学者会直接拿二阶方程去套结果报错。这里的核心步骤是用状态变量替换令x1xx2x_dot那么原方程就能改写成$$\begin{cases} \frac{dx_1}{dt}x_2 \ \frac{dx_2}{dt}\frac{F(t)-cx_2-kx_1}{m} \end{cases}$$这种写法在MATLAB中是一个长度为2的列向量dx。ODE45只负责在自适应步长上计算x1和x2的数值至于怎么算dx完全由你提供的函数决定。所以理解这个状态空间转换是完成求解的前提。我见过不少人卡在这一步其实只要把“位移”和“速度”当作两个独立变量问题就迎刃而解。2.2 MATLAB函数编写与ODE45调用下面是一个可以直接改的MATLAB函数示例。我把模型参数和故障力函数分开写方便调试function dx bearing_ode(t, x, model, fault_force) x1 x(1); x2 x(2); m model.m; c model.c; k model.k; F fault_force(t); dx zeros(2,1); dx(1) x2; dx(2) (F - c*x2 - k*x1) / m; end故障力函数里我预先计算好故障冲击时刻序列然后在每个冲击时刻邻域内生成衰减振荡。一个简单的实现是function F fault_force(t) global fault_times A zeta_f wn_f F 0; for i 1:length(fault_times) dt t - fault_times(i); if dt 0 dt 0.01 F F A * exp(-zeta_f*wn_f*dt) * sin(wn_f*sqrt(1-zeta_f^2)*dt); end end end当然用global变量不是好习惯实际项目中建议用嵌套函数或者参数对象传递。但核心思路不变。调用ODE45时代码如下tspan [0 1]; % 仿真1秒 x0 [0; 0]; % 初始位移和速度 options odeset(RelTol,1e-6,AbsTol,1e-8); [t, x] ode45((t,x) bearing_ode(t,x,model,fault_force), tspan, x0, options);这里options设置很重要。默认容差下ODE45可能为了满足精度要求在冲击附近疯狂缩小步长导致计算量暴增。我通常把RelTol设在1e-6左右AbsTol设在1e-8左右既能保证精度又不至于太慢。2.3 ODE45误差控制与信号重采样ODE45是变步长求解器输出时间点不是均匀的。可我们在后续做FFT和包络谱时需要等间隔采样的信号。所以求解完成后必须用interp1把结果重采样到目标采样率fs上。例如把输出插值到fs20000 Hz对应的均匀时间网格上fs 20000; t_uniform (0:1/fs:tspan(end)); x_uniform interp1(t, x(:,1), t_uniform, spline);这里插值方法选用‘spline’是因为三次样条在光滑性上比线性插值好不会引入额外的高频伪迹。如果你的仿真时长是1秒采样率2万赫兹那就会得到2万个均匀采样点足够覆盖故障特征频率的几十次谐波。我实测下来当模型刚度和冲击频率都处在合理范围时ODE45在1秒仿真时长内大概能在几秒内完成。如果发现计算慢优先检查是不是冲击间隔太短导致ODE45在极短时间步长里反复计算。这时候可以把冲击力宽度适当加大或者转用刚性求解器后面会详细讲。3. 故障特征怎么提取时域、频域与包络谱3.1 仿真信号的故障特征频率计算要验证仿真信号对不对最关键的是计算故障特征频率然后看包络谱峰值是否落在这些频率上。滚动轴承的典型故障特征频率公式如下外圈故障特征频率$$BPFO \frac{N_b f_r}{2}\left(1 - \frac{d}{D}\cos\alpha\right)$$内圈故障特征频率$$BPFI \frac{N_b f_r}{2}\left(1 \frac{d}{D}\cos\alpha\right)$$滚动体故障特征频率$$BSF \frac{D f_r}{2d}\left[1 - \left(\frac{d}{D}\cos\alpha\right)^2\right]$$其中N_b是滚动体个数f_r是转轴频率d是滚动体直径D是节径alpha是接触角。以我常用的示例参数为例滚动体数N_b9转频f_r30Hz节径D65mm滚动体直径d15mm接触角alpha0度。代入公式Nb 9; d 15e-3; D 65e-3; alpha 0; fr 30; BPFO (Nb*fr/2)*(1 - d/D*cos(alpha)); BPFI (Nb*fr/2)*(1 d/D*cos(alpha)); BSF (D*fr/(2*d))*(1 - (d/D*cos(alpha))^2); disp([BPFO, BPFI, BSF]);运行后得到外圈故障频率约103.85Hz内圈故障频率约166.15Hz滚动体故障频率约81.05Hz。这个数值在后面的包络谱上应该是非常明显的峰值。仿真时冲击间隔就按103.85Hz对应的周期产生也就是每隔约9.6毫秒产生一次冲击。3.2 从时域冲击到包络谱的具体实现先看时域波形。把ODE45求出的位移做一次差分可以得到加速度波形或者在ODE函数里直接以加速度作为第三个输出项。时域波形中外圈故障会有明显的等间隔冲击内圈故障的冲击幅值会随着承载区变化而调制这也是区分故障类型的一个重要线索。然后做包络谱。流程是带通滤波-希尔伯特变换取包络-对包络信号做FFT。带通滤波的目的是去掉低频转频成分和高频噪声只保留系统共振频带。这里带通中心频率可以设置成与模型的固有频率一致。例如我用2000Hz到10000Hz的带通区间。% 假设信号是均匀重采样后的加速度accfs20000 [b, a] butter(2, [2000 10000]/(fs/2), bandpass); acc_filt filtfilt(b, a, acc); env abs(hilbert(acc_filt)); % 希尔伯特包络 env env - mean(env); % 去掉直流分量 N length(env); f_axis (0:N-1) * (fs/N); env_spectrum abs(fft(env));绘制包络谱时重点关注BPFO及其2倍频、3倍频。如果在这些位置出现明显峰值可以判断为外圈故障。同理用内圈故障冲击序列重跑仿真看BPFI峰值即可。这样仿真数据和诊断算法就形成了一个闭环哪一步出了问题都能定位。3.3 仿真到诊断的流程闭环这个流程最大的优势是“可重复可控制”。你可以通过改变冲击幅值A、冲击衰减系数zeta_f、故障特征频率等参数模拟不同故障程度和不同工况。比如外圈故障早期冲击幅值小、噪声大严重故障时冲击间隔可能变得不稳定。你可以把这种变化注入到仿真信号里用来测试你的诊断算法在低信噪比下是否依然有效。我自己的做法是先把外圈、内圈、滚动体三种故障各生成一批仿真数据然后跑一个简单的自动诊断脚本计算包络谱再提取BPFO/BPFI/BSF附近的谱峰幅度超过阈值就判为对应故障。这套逻辑虽然在真实数据上还需要再校准阈值但算法框架的验证已经足够了。后续换真实数据时需要调整的只是滤波参数和阈值骨干逻辑不用改。4. 常见问题与排查技巧实录4.1 ODE45刚性问题与刚性求解器很多人在跑轴承模型时发现ODE45运行时间特别长甚至卡死。这种情况大概率是因为模型出现了“刚性问题”。所谓刚性通俗讲就是解中包含的快变分量和慢变分量差别太大。比如接触刚度k取了非常大的值系统固有频率达到十几千赫兹而故障冲击间隔长达几十毫秒ODE45为了保证快变的精度不得不把步长压得极小自然就慢。解决方式有两种。最直接的是换用ode15s或ode23t这类刚性求解器。虽然在标准轴承模型里非刚性算法通常也能跑但在参数极端时刚性求解器会快很多。另一种方式是从建模角度调整参数例如适当减小刚度使系统固有频率靠近实际传感器可测范围几千赫兹这样ODE45就不会因为高频分量耗费过多时间。如果两种方法都不适用还可以对时间变量做无量纲化处理把快慢分量的尺度拉近。4.2 仿真结果发散或震荡有时仿真结果在某个时刻突然变成NaN或者指数增长原因不外乎初始位移或速度过大、激励力幅值太大、阻尼比设成了负数。这里有个经验小贴士先不加载故障激励只给系统一个初始位移看它是否能自然衰减到零。如果衰减正常说明阻尼和刚度设置没问题如果发散优先检查参数是否符号错误。然后再逐步加大故障力幅值看响应是否可控。这种“分步排查”法比对着错误信息瞎猜快得多。另外一个容易被忽略的问题是ODE45在冲击力函数的间断点上可能“感受不到”冲击。如果fault_force只在某个时刻返回非零值而ODE45恰好在那个时刻没有取点冲击就会被漏掉。解决办法是设置MaxStep把最大步长限制在冲击时长的1/10以内或者给冲击力加一个小斜坡让它不要那么突变。我一般会给冲击力加一个很短的上升沿既符合物理过程也能避免数值问题。4.3 包络谱看不出故障频率怎么办包络谱里没看到BPFO峰值先不急着改算法按顺序排查先看时域信号是否存在周期性冲击然后用cpsd或者简单FFT看共振频带能量是否明显再看带通滤波频率区间是否对准了系统固有频率。很多初学者把带通区间设置得很随意结果把故障信息全滤掉了。还有一个常见坑是采样率不足。如果仿真时重采样率fs太低计算出的FFT频率分辨率不够或者奈奎斯特频率太小故障特征频率的高次谐波会被混叠。建议fs至少设为最大关注频率的5到10倍。例如如果你希望在10kHz以内分析fs至少要50kHz我通常取20kHz能覆盖最高共振频率约9kHz效果已经不错。最后检查一下冲击序列是否确实存在。用findpeaks函数在包络信号里找冲击间隔对比理论故障周期这个动作看起来简单但能在5秒内定位问题。下面用一个小表格总结常见现象和排查方向现象可能原因排查顺序ODE45运行极慢模型刚性/冲击步长太小换ode15s降刚度设MaxStep结果发散/NaN阻尼或刚度参数异常取消激励看自由响应逐步加大激励时域没有冲击故障力函数未触发检查冲击序列和MaxStep包络谱无峰带通滤波区间偏移先看FFT共振频带再调整滤波器峰位置偏差特征频率公式参数错误重新核对轴承几何参数和f_r我在实际使用中还有一个习惯每次仿真都会打印BPFO、BPFI这些理论值然后在包络谱图上用垂直线标记出来。这样即使诊断算法出问题也能第一时间判断是仿真数据的问题还是特征提取流程的问题。调试效率提升非常明显。最后再分享一个小技巧。如果你初期只是想验证ODE45和包络谱流程不要一上来就追求复杂的非线性接触模型。先用最简单的线性质量-弹簧-阻尼系统把外圈故障跑通再逐步加入内圈、滚动体故障甚至加入随机噪声和转速波动。这套“由简到繁”的玩法能让每个环节的坑都暴露得明明白白也方便你对照理论公式验证代码正确性。我现在每次接到新的诊断需求都会先回到这个基础仿真框架再根据实际信号特征迭代模型复杂度这比直接上手写一堆复杂代码要稳得多。本文还有配套的精品资源点击获取 
长河君 · 2026/8/31 8:05:27 · 阅读 1.2k MATLAB+ODE45实现轴承故障动力学建模与包络谱诊断
简介本资源面向机械工程、故障诊断与振动分析领域的高校师生及工业设备维护工程师聚焦轴承动力学建模与早期故障识别这一核心问题提供基于MATLAB的可运行仿真与诊断方案。压缩包共6个文件全部为.m脚本含vdp类与v类命名模型文件总大小仅3KB轻量紧凑涵盖轴承非线性动力学方程构建、ODE45高精度数值求解、振动响应仿真及典型故障特征提取等关键环节代码结构清晰、注释友好便于理解滚动体-套圈接触力学建模逻辑与故障信号演化规律。已有3009人学习下载适用于课程设计、科研建模入门及预测性维护算法验证场景读者可直接调用脚本复现轴承在正常、内圈/外圈缺陷等工况下的时频域动态响应并结合信号处理思路开展峭度分析、谱图判读等诊断实践是贯通理论建模与工程诊断的实用型MATLAB教学与开发参考。 做轴承故障诊断的工程人员和研究生基本都逃不过一个坎手里只有一堆振动数据想验证诊断算法却拿不到带故障的轴承样本。买一套故障轴承不便宜试验台搭建周期又长而且故障程度还不一定可控。我的做法是用MATLAB先做轴承动力学建模用ODE45把运动方程解出来仿真出一组“带故障”的振动信号再拿它调试包络谱、特征频率识别这些诊断逻辑。这套流程我用了很多次稳定、可控而且能重复生成任意工况的数据比纯理论推导直观也比纯数据驱动的方法更容易定位问题。这篇文章就完整记录一下从模型建立、参数设置、ODE45求解到故障诊断代码实现希望给正在做相关方向的朋友提供一个可以直接参考的起点。1. 轴承动力学模型从哪来问题定义与建模思路1.1 轴承故障诊断为什么需要动力学模型直接回答一个问题既然有振动信号为什么不直接做频谱分析非要先建模我的理解是实际采集到的轴承信号里混着转频、齿轮啮合频率、环境噪声甚至传感器安装共振很难干净地把故障特征分离出来。动力学模型的价值在于它把“故障”变成一个可量化的激励源你可以控制故障位置、故障程度、转速和负载然后观察振动响应的变化规律。这个过程中建立的“特征频率-故障类型”对应关系比单纯从数据里翻出来的规律可靠得多也更容易解释。另外在开发诊断算法时真实故障数据是稀缺资源。你不可能为了调一个滤波器参数就去买几十个故障轴承。仿真模型相当于一个“万能试验台”你可以批量生成数据然后测试不同的诊断方法时域指标、包络谱、倒频谱、深度学习分类器都能在这套数据上先跑通。我自己的习惯是先用仿真数据验证算法逻辑再拿实际数据测试这样能省下大量调试时间。1.2 单自由度质量-弹簧-阻尼模型够不够用最早我尝试过直接用有限元模型做一个完整的三维滚动轴承网格一划计算一次要好几个小时而且参数调整极不灵活。后来慢慢发现对于故障诊断这种“看特征频率”的场景单自由度模型已经完全够用。它的物理意义非常清晰把轴承-转子系统简化为一个等效质量m通过一个等效刚度k和阻尼c支承在基础上。当滚动体经过局部缺陷时会产生一个周期性冲击力F(t)系统的振动响应就是对这个冲击的衰减振荡。最基本的运动方程写出来是$$m\ddot{x}c\dot{x}kxF(t)$$其中x可以理解为轴承外圈或传感器安装位置的等效位移。m、c、k不是直接从轴承手册里查到的而是根据系统共振行为反推的等效值。实际操作中先估计系统的固有频率fn和阻尼比zeta然后令km*(2pifn)^2c2zetasqrt(k*m)。这样做的好处是你只需要关注你关心的频带模型动态特性基本可控。对于包络谱分析我们通常关注的是冲击激起的系统高频共振所以固有频率可以设置在1kHz到5kHz之间这也是很多实际诊断中带通滤波器常用的频带。1.3 故障冲击激励怎么模拟才真实故障激励不是纯粹的脉冲函数因为滚动体通过缺陷时接触状态是渐变的。常用的做法是用一个带高频衰减振荡的短时冲击来表示。也就是说每次冲击不是理想冲激而是激发系统固有频率后快速衰减的振荡。如果直接用一个Dirac函数ODE45会因为阶跃过于尖锐而需要极小的步长算得很慢也不符合物理过程。我通常把故障力建模成如下形式$$F_f(t)\sum_{i} A_i \exp\left[-\zeta_f \omega_n (t-t_i)\right]\sin\left(\omega_n \sqrt{1-\zeta_f^2}(t-t_i)\right) \cdot u(t-t_i)$$其中t_i是第i次冲击的发生时刻A_i是冲击幅值zeta_f控制冲击衰减速度omega_n是冲击响应的角频率。这个形式其实模拟了“滚动体撞击缺陷边缘-激发局部共振-逐渐衰减”的过程。冲击时刻t_i按照故障特征频率的倒数间隔产生也就是每经过一个缺陷产生一次冲击。为了让仿真更贴近真实信号我还会在激励上加入小幅度的随机扰动模拟转速波动和载荷变化。2. 用ODE45求解运动方程从数学方程到数值解2.1 二阶方程转一阶状态空间ODE45是MATLAB中基于Dormand-Prince方法的非刚性常微分方程求解器它要求你把待求解的方程写成一阶微分方程组的形式。很多初学者会直接拿二阶方程去套结果报错。这里的核心步骤是用状态变量替换令x1xx2x_dot那么原方程就能改写成$$\begin{cases} \frac{dx_1}{dt}x_2 \ \frac{dx_2}{dt}\frac{F(t)-cx_2-kx_1}{m} \end{cases}$$这种写法在MATLAB中是一个长度为2的列向量dx。ODE45只负责在自适应步长上计算x1和x2的数值至于怎么算dx完全由你提供的函数决定。所以理解这个状态空间转换是完成求解的前提。我见过不少人卡在这一步其实只要把“位移”和“速度”当作两个独立变量问题就迎刃而解。2.2 MATLAB函数编写与ODE45调用下面是一个可以直接改的MATLAB函数示例。我把模型参数和故障力函数分开写方便调试function dx bearing_ode(t, x, model, fault_force) x1 x(1); x2 x(2); m model.m; c model.c; k model.k; F fault_force(t); dx zeros(2,1); dx(1) x2; dx(2) (F - c*x2 - k*x1) / m; end故障力函数里我预先计算好故障冲击时刻序列然后在每个冲击时刻邻域内生成衰减振荡。一个简单的实现是function F fault_force(t) global fault_times A zeta_f wn_f F 0; for i 1:length(fault_times) dt t - fault_times(i); if dt 0 dt 0.01 F F A * exp(-zeta_f*wn_f*dt) * sin(wn_f*sqrt(1-zeta_f^2)*dt); end end end当然用global变量不是好习惯实际项目中建议用嵌套函数或者参数对象传递。但核心思路不变。调用ODE45时代码如下tspan [0 1]; % 仿真1秒 x0 [0; 0]; % 初始位移和速度 options odeset(RelTol,1e-6,AbsTol,1e-8); [t, x] ode45((t,x) bearing_ode(t,x,model,fault_force), tspan, x0, options);这里options设置很重要。默认容差下ODE45可能为了满足精度要求在冲击附近疯狂缩小步长导致计算量暴增。我通常把RelTol设在1e-6左右AbsTol设在1e-8左右既能保证精度又不至于太慢。2.3 ODE45误差控制与信号重采样ODE45是变步长求解器输出时间点不是均匀的。可我们在后续做FFT和包络谱时需要等间隔采样的信号。所以求解完成后必须用interp1把结果重采样到目标采样率fs上。例如把输出插值到fs20000 Hz对应的均匀时间网格上fs 20000; t_uniform (0:1/fs:tspan(end)); x_uniform interp1(t, x(:,1), t_uniform, spline);这里插值方法选用‘spline’是因为三次样条在光滑性上比线性插值好不会引入额外的高频伪迹。如果你的仿真时长是1秒采样率2万赫兹那就会得到2万个均匀采样点足够覆盖故障特征频率的几十次谐波。我实测下来当模型刚度和冲击频率都处在合理范围时ODE45在1秒仿真时长内大概能在几秒内完成。如果发现计算慢优先检查是不是冲击间隔太短导致ODE45在极短时间步长里反复计算。这时候可以把冲击力宽度适当加大或者转用刚性求解器后面会详细讲。3. 故障特征怎么提取时域、频域与包络谱3.1 仿真信号的故障特征频率计算要验证仿真信号对不对最关键的是计算故障特征频率然后看包络谱峰值是否落在这些频率上。滚动轴承的典型故障特征频率公式如下外圈故障特征频率$$BPFO \frac{N_b f_r}{2}\left(1 - \frac{d}{D}\cos\alpha\right)$$内圈故障特征频率$$BPFI \frac{N_b f_r}{2}\left(1 \frac{d}{D}\cos\alpha\right)$$滚动体故障特征频率$$BSF \frac{D f_r}{2d}\left[1 - \left(\frac{d}{D}\cos\alpha\right)^2\right]$$其中N_b是滚动体个数f_r是转轴频率d是滚动体直径D是节径alpha是接触角。以我常用的示例参数为例滚动体数N_b9转频f_r30Hz节径D65mm滚动体直径d15mm接触角alpha0度。代入公式Nb 9; d 15e-3; D 65e-3; alpha 0; fr 30; BPFO (Nb*fr/2)*(1 - d/D*cos(alpha)); BPFI (Nb*fr/2)*(1 d/D*cos(alpha)); BSF (D*fr/(2*d))*(1 - (d/D*cos(alpha))^2); disp([BPFO, BPFI, BSF]);运行后得到外圈故障频率约103.85Hz内圈故障频率约166.15Hz滚动体故障频率约81.05Hz。这个数值在后面的包络谱上应该是非常明显的峰值。仿真时冲击间隔就按103.85Hz对应的周期产生也就是每隔约9.6毫秒产生一次冲击。3.2 从时域冲击到包络谱的具体实现先看时域波形。把ODE45求出的位移做一次差分可以得到加速度波形或者在ODE函数里直接以加速度作为第三个输出项。时域波形中外圈故障会有明显的等间隔冲击内圈故障的冲击幅值会随着承载区变化而调制这也是区分故障类型的一个重要线索。然后做包络谱。流程是带通滤波-希尔伯特变换取包络-对包络信号做FFT。带通滤波的目的是去掉低频转频成分和高频噪声只保留系统共振频带。这里带通中心频率可以设置成与模型的固有频率一致。例如我用2000Hz到10000Hz的带通区间。% 假设信号是均匀重采样后的加速度accfs20000 [b, a] butter(2, [2000 10000]/(fs/2), bandpass); acc_filt filtfilt(b, a, acc); env abs(hilbert(acc_filt)); % 希尔伯特包络 env env - mean(env); % 去掉直流分量 N length(env); f_axis (0:N-1) * (fs/N); env_spectrum abs(fft(env));绘制包络谱时重点关注BPFO及其2倍频、3倍频。如果在这些位置出现明显峰值可以判断为外圈故障。同理用内圈故障冲击序列重跑仿真看BPFI峰值即可。这样仿真数据和诊断算法就形成了一个闭环哪一步出了问题都能定位。3.3 仿真到诊断的流程闭环这个流程最大的优势是“可重复可控制”。你可以通过改变冲击幅值A、冲击衰减系数zeta_f、故障特征频率等参数模拟不同故障程度和不同工况。比如外圈故障早期冲击幅值小、噪声大严重故障时冲击间隔可能变得不稳定。你可以把这种变化注入到仿真信号里用来测试你的诊断算法在低信噪比下是否依然有效。我自己的做法是先把外圈、内圈、滚动体三种故障各生成一批仿真数据然后跑一个简单的自动诊断脚本计算包络谱再提取BPFO/BPFI/BSF附近的谱峰幅度超过阈值就判为对应故障。这套逻辑虽然在真实数据上还需要再校准阈值但算法框架的验证已经足够了。后续换真实数据时需要调整的只是滤波参数和阈值骨干逻辑不用改。4. 常见问题与排查技巧实录4.1 ODE45刚性问题与刚性求解器很多人在跑轴承模型时发现ODE45运行时间特别长甚至卡死。这种情况大概率是因为模型出现了“刚性问题”。所谓刚性通俗讲就是解中包含的快变分量和慢变分量差别太大。比如接触刚度k取了非常大的值系统固有频率达到十几千赫兹而故障冲击间隔长达几十毫秒ODE45为了保证快变的精度不得不把步长压得极小自然就慢。解决方式有两种。最直接的是换用ode15s或ode23t这类刚性求解器。虽然在标准轴承模型里非刚性算法通常也能跑但在参数极端时刚性求解器会快很多。另一种方式是从建模角度调整参数例如适当减小刚度使系统固有频率靠近实际传感器可测范围几千赫兹这样ODE45就不会因为高频分量耗费过多时间。如果两种方法都不适用还可以对时间变量做无量纲化处理把快慢分量的尺度拉近。4.2 仿真结果发散或震荡有时仿真结果在某个时刻突然变成NaN或者指数增长原因不外乎初始位移或速度过大、激励力幅值太大、阻尼比设成了负数。这里有个经验小贴士先不加载故障激励只给系统一个初始位移看它是否能自然衰减到零。如果衰减正常说明阻尼和刚度设置没问题如果发散优先检查参数是否符号错误。然后再逐步加大故障力幅值看响应是否可控。这种“分步排查”法比对着错误信息瞎猜快得多。另外一个容易被忽略的问题是ODE45在冲击力函数的间断点上可能“感受不到”冲击。如果fault_force只在某个时刻返回非零值而ODE45恰好在那个时刻没有取点冲击就会被漏掉。解决办法是设置MaxStep把最大步长限制在冲击时长的1/10以内或者给冲击力加一个小斜坡让它不要那么突变。我一般会给冲击力加一个很短的上升沿既符合物理过程也能避免数值问题。4.3 包络谱看不出故障频率怎么办包络谱里没看到BPFO峰值先不急着改算法按顺序排查先看时域信号是否存在周期性冲击然后用cpsd或者简单FFT看共振频带能量是否明显再看带通滤波频率区间是否对准了系统固有频率。很多初学者把带通区间设置得很随意结果把故障信息全滤掉了。还有一个常见坑是采样率不足。如果仿真时重采样率fs太低计算出的FFT频率分辨率不够或者奈奎斯特频率太小故障特征频率的高次谐波会被混叠。建议fs至少设为最大关注频率的5到10倍。例如如果你希望在10kHz以内分析fs至少要50kHz我通常取20kHz能覆盖最高共振频率约9kHz效果已经不错。最后检查一下冲击序列是否确实存在。用findpeaks函数在包络信号里找冲击间隔对比理论故障周期这个动作看起来简单但能在5秒内定位问题。下面用一个小表格总结常见现象和排查方向现象可能原因排查顺序ODE45运行极慢模型刚性/冲击步长太小换ode15s降刚度设MaxStep结果发散/NaN阻尼或刚度参数异常取消激励看自由响应逐步加大激励时域没有冲击故障力函数未触发检查冲击序列和MaxStep包络谱无峰带通滤波区间偏移先看FFT共振频带再调整滤波器峰位置偏差特征频率公式参数错误重新核对轴承几何参数和f_r我在实际使用中还有一个习惯每次仿真都会打印BPFO、BPFI这些理论值然后在包络谱图上用垂直线标记出来。这样即使诊断算法出问题也能第一时间判断是仿真数据的问题还是特征提取流程的问题。调试效率提升非常明显。最后再分享一个小技巧。如果你初期只是想验证ODE45和包络谱流程不要一上来就追求复杂的非线性接触模型。先用最简单的线性质量-弹簧-阻尼系统把外圈故障跑通再逐步加入内圈、滚动体故障甚至加入随机噪声和转速波动。这套“由简到繁”的玩法能让每个环节的坑都暴露得明明白白也方便你对照理论公式验证代码正确性。我现在每次接到新的诊断需求都会先回到这个基础仿真框架再根据实际信号特征迭代模型复杂度这比直接上手写一堆复杂代码要稳得多。本文还有配套的精品资源点击获取