干设备振动监测这一行打交道最多的就是滚动轴承尤其是深沟球轴承。早期我在某实验室做故障诊断课题时被一个很现实的问题卡住了手上有一堆现场采集的轴承振动数据想验证自己设计的诊断算法但正常理论解析又没有解析解拿实验数据一对照又总被各种工况噪声干扰。后来一个偶然机会接触到Matlab仿真建模发现用集中质量模型加Hertz接触理论居然能在个人电脑上复现轴承振动响应的核心特征而且成本极低。这正是这篇内容的价值所在不需要昂贵的实验台不依赖复杂的有限元软件只要会基本的Matlab语法就能把滚动轴承的时变刚度、局部缺陷激励、非线性接触这些核心动力学行为跑出来。这篇内容适合谁如果你是刚接触转子动力学和故障诊断的研究生或者在公司做设备状态监测的工程师想搞懂特征频率BPFO、BPFI到底是怎么在振动信号里“长”出来的这篇文章能帮你把理论公式和仿真代码一一对上。它不追求像有限元那样对每个接触节点的应力都精细刻画而是聚焦在宏观动力学响应——即轴承内圈相对外圈的振动位移轨迹这正是事后包络分析、频域特征提取的基础。你跟着操作一遍不仅能拿到能跑的Matlab代码还能理解背后的参数选择逻辑和常见的坑。1. 建模前的思路先想清楚仿真对象和输出量1.1 这个模型到底在算什么深沟球轴承的动力学建模往细了说可以做得极其复杂——滚动体与套圈的弹性流体动压润滑、保持架与滚动体之间的碰撞、热变形耦合等等。但对初学者最重要的是建立一个能抓住主要矛盾的低自由度模型。通常的做法是“集中质量法”把内圈和转轴的组合质量等效为一个集中质量块外圈固定在轴承座或壳体上视为刚体滚动体被当作只有刚度效应、没有质量的非线性弹簧。那么我们要跟踪的状态量就很简单内圈质量块在水平方向X方向和垂直方向Y方向的位移及速度。为什么要这么简化因为对故障诊断来说我们关心的不是轴承内部某个点的应力而是内外圈之间的相对运动会如何通过轴承座传递到壳体上形成传感器能测到的振动加速度。自由度越多参数标定越痛苦而顶部2自由度模型已经能捕捉滚动体通过频率及其谐波足够解释频谱上大部分谱峰来源。往外再扩展一下如果关心倾斜力矩的影响可以再加一个绕X轴或Y轴的旋转自由度如果考虑转轴弯曲可以做成多节点转子轴承系统。但这一步一步来先把2自由度模型跑通、跑出正确频率特征后面做扩展就顺理成章。1.2 定下驱动这场振动的“剧本”转频与特征频率轴承振动不是随机噪声滚动体在滚道内滚动时会以固定的规律周期性通过某一点。这个“周期”用频率来表达就是特征频率。建模前算清楚这些频率不仅是为了后面验证频谱峰值更是为了在设计仿真时长和采样率时心里有数。深沟球轴承的主要特征频率如下接触角近似为0°外圈故障频率BPFO( BPFO \frac{n \cdot f_r}{2} \left(1 - \frac{d_b}{D} \right) )内圈故障频率BPFI( BPFI \frac{n \cdot f_r}{2} \left(1 \frac{d_b}{D} \right) )滚动体故障频率BSF( BSF \frac{D \cdot f_r}{2 \cdot d_b} \left(1 - \left(\frac{d_b}{D}\right)^2 \right) )保持架故障频率FTF( FTF \frac{f_r}{2} \left(1 - \frac{d_b}{D} \right) )其中 ( n ) 是滚动体数量( f_r ) 是转轴转频Hz( d_b ) 是滚动体直径( D ) 是节圆直径pitch diameter。现场做诊断时我们用这些公式判断频谱上的峰属于哪一类缺陷做仿真时我们会人为在某个位置设置缺陷确保仿真结果在对应频率处出现峰值这就完成了理论到信号的闭环。2. 核心原理拆解Hertz接触、时变刚度和间隙2.1 接触刚度的非线性球轴承滚动体与滚道之间的接触属于典型的点接触力学上通常用Hertz接触理论描述。滚动体压入滚道时的法向载荷 ( Q ) 与法向变形量 ( \delta ) 之间满足[ Q K \delta^{3/2} ]注意这里的指数3/2它来自Hertz点接触的椭圆接触面积与载荷的非线性关系。这个非线性是轴承振动区别于普通线性弹簧系统的核心原因。你如果把它简化为线性弹簧 ( Q K\delta )那么仿真谱上很多高次谐波和幅值跳跃现象会完全消失。那么接触刚度系数 ( K ) 怎么来严格计算出它需要通过滚动体和套圈的曲率半径、材料弹性模量和泊松比查表或迭代。对于初学者我建议先用一个经验量级钢制深沟球轴承( K ) 量级通常在 ( 10^8 \sim 10^9 , \text{N/m}^{3/2} ) 之间。启动时可以取 ( K 1 \times 10^8 )跑出结果后再做参数敏感性分析。这比一上来就查表展开椭圆积分要友好得多。实际写代码时还有个细节( \delta ) 是滚动体与滚道的总变形量不是说位移代进去就是变形。接触变形 ( \delta_j ) 需要根据该滚动体位置的几何关系来算当该位置滚动体受压时变形量为正若间隙使滚动体脱离接触变形量会被截断为零。这就是下一节要讲的关键几何。2.2 滚动体位置与系统周期性的来源深沟球轴承转动时滚动体绕轴公转。设保持架转速为 ( \omega_c 2\pi \cdot FTF )那么第 ( j ) 个滚动体的角位置可以写成[ \theta_j \omega_c t \theta_0 \frac{2\pi (j-1)}{n} ]为什么要关心角位置因为接触力分解到X、Y方向时需要知道每个滚动体的受力方向。而且当滚动体个数为 ( n ) 时系统对外呈现的“通过特征”是 ( n ) 倍保持架频率——在每个滚动体经过承载区时支撑刚度都会发生一次周期波动这就是时变刚度的源头。有一个很容易踩的坑初始时刻的角位置 ( \theta_0 ) 不是随便给的它会直接影响仿真的瞬态过程和某些谐波相位。建议在有实验基准时把 ( \theta_0 ) 与位置传感器的初始测量同步没有基准时随意给一个角度也能跑只是瞬态段长一点。2.3 径向间隙和润滑阻尼怎么取舍实际轴承内部有径向游隙承载区只在下方若干滚动体之间分布。仿真中这个“游隙”非常关键因为它是非线性重力和接触力的“开关”只有当总位移让某个滚动体受压时该滚动体才贡献刚度。这个开关会导致滚动体进出承载区时刚度猛变产生冲击性振动。模型中把初始径向间隙记为 ( c_0 )单位米。考虑重力方向通常沿负Y方向在内圈自重的压缩下内圈会下沉一段距离底部几个滚动体先接触其他滚动体与滚道之间保留间隙。仿真时每个滚动体的总变形量[ \delta_j x \cos \theta_j y \sin \theta_j - c_0 ]注意符号约定( x, y ) 是内圈质心的位移。但要小心不同文献坐标系定义不同有的把 ( \theta_j ) 从X轴起算有的从负Y轴起算。我在自己调试时踩过这个坑符号反了的结果是振动激励方向全反频谱照样出现特征频率但时域波形会差一个相位初学阶段很难自我发现。稳妥做法是先做一组无缺陷仿真用解析解或文献结果对比验证坐标方向。阻尼方面滚动体-滚道接触处的阻尼机制非常复杂。工程上通常用等效阻尼比近似在内圈运动方程里加一个线性阻尼项例如取阻尼比 ( \zeta 0.01 \sim 0.05 )。模型代码中对应的阻尼系数是 ( C 2\zeta \sqrt{K_dyn m} )这里的 ( K_dyn ) 可以用平均接触刚度的线性化值不必精确。3. 逐步搭建Matlab模型从参数到求解3.1 参数定义与初始化我习惯把所有参数集中放在一个脚本开头方便改工况。先给一台某实验室用的单支承轴承做仿真参数如下滚动体数量 ( n 9 )滚动体直径 ( d_b 7.94, \text{mm} )节圆直径 ( D 38.5, \text{mm} )径向间隙 ( c_0 20, \mu m )转频 ( f_r 30, \text{Hz} )对应1800 rpm内圈轴等效质量 ( m 1.2, \text{kg} )接触刚度 ( K 1 \times 10^8 )阻尼比 ( \zeta 0.02 )Matlab里注意单位全部换算成国际单位米、千克、秒这个不能偷懒。好多人调试出天文数字的加速度就是毫米和米混用。初始化阶段除了定义参数还要设置仿真时间和采样输出需求。ode45会自动变步长但后处理做FFT时需要等间隔数据所以仿真完要用interp1重采样。不妨一开始就规划好仿真时长0.5秒重采样频率为10 kHz覆盖到几百Hz的特征频率和调制边带足够。3.2 组装运动方程与接触力函数状态向量设为 ( [x, y, v_x, v_y] )那么加速度[ m \ddot{x} c \dot{x} F_{cx} F_{u,x} ] [ m \ddot{y} c \dot{y} F_{cy} - mg F_{u,y} ]这里的 ( F_{cx}, F_{cy} ) 是所有滚动体接触力在X、Y方向上的合力方向取决于滚动体位置。( F_{u,x}, F_{u,y} ) 是外部激励。无外部载荷时( F_{u,x}0 )( F_{u,y}0 )但自重项 (-mg) 不能丢因为它是让轴承进入承载区的“启动动力”。求解接触力时要遍历每个滚动体先算角位置再算变形量截负为0然后[ F_j K \cdot \delta_j^{3/2} ] [ F_{cx} -\sum_j F_j \cos \theta_j ] [ F_{cy} -\sum_j F_j \sin \theta_j ]为什么合力是负号因为接触力是把滚动体压缩的内圈往滚动体方向拉等效为抵抗内圈位移的恢复力。用Matlab写一个嵌套函数或独立函数function dstate bearingDynamics(t, state, params) x state(1); y state(2); vx state(3); vy state(4); Fcx 0; Fcy 0; for j 1:params.n theta_j params.wc * t params.theta0 2*pi*(j-1)/params.n; delta x * cos(theta_j) y * sin(theta_j) - params.c0; if delta 0 F params.K * delta ^ 1.5; Fcx Fcx - F * cos(theta_j); Fcy Fcy - F * sin(theta_j); end end dstate zeros(4,1); dstate(1) vx; dstate(2) vy; dstate(3) (Fcx params.Fux) / params.m; dstate(4) (Fcy - params.m * params.g params.Fuy) / params.m; end这段代码是核心初学者不要照抄后就不管了。你要能自己手工推导一次合力表达式确认坐标的正负号和你预设的滚动体旋转方向一致。3.3 设置初始位移并调用ode45这是个看起来小但实际很关键的环节。如果把初始位移设成0那么所有滚动体初始都处于间隙状态重力作用下内圈自由落体接触力从零猛增仿真初期的瞬态剧烈振荡会拉长稳定时间有时甚至导致ode45报错。稳妥的做法是预估一个静态下沉量。简化计算假设所有底部滚动体均匀承载总的载荷 ( m g F_{load} ) 由 ( n/2 ) 个滚动体分担那么每个滚动体的变形可以通过 ( K \delta^{3/2} \text{载荷}/n_{eff} ) 反算。例如总载荷约11.8 N有效滚动体3个每个约4 N( \delta \approx (4/1e8)^{2/3} \approx 1.1 \times 10^{-5} , \text{m} )也就是11微米。初始位移设置为 ( x0, y-\delta_{est} )向重力方向压缩能明显减少瞬态。只要模型是稳定的初始偏移偏离真实值也没关系前0.05秒瞬态会被阻尼吸收取后0.4秒做分析即可。调用求解器我可以直接给出示例params.g 9.81; params.wc 2*pi* (30/2) * (1 - 7.94e-3/38.5e-3); params.theta0 0; tspan [0 0.5]; state0 [0; -1.1e-5; 0; 0]; [t, state] ode45((t, s) bearingDynamics(t, s, params), tspan, state0);需要注意区分轴转频 ( f_r ) 和保持架转速 ( f_{FTF} )( \omega_c ) 用的是 ( 2\pi FTF )不要拿轴转频直接算滚动体位置否则时变刚度的特征频率会整体偏高。求解完成后用线性插值重采样fs 10000; tq linspace(0, 0.5, 0.5*fs 1); Xq interp1(t, state(:, 1), tq); Yq interp1(t, state(:, 2), tq);取后半段做频谱分析就能看到转频和滚动体通过频率的峰值。如果时域波形是杂乱无章的先不要急着调参数检查两件事一是滚动体总个数及初始角位置是否导致系统周期性被破坏二是阻尼是否太小导致在0.5秒内尚未收敛到稳态。4. 进阶玩法给轴承“制造”一个局部缺陷4.1 缺陷建模的常见套路健康轴承的动力学仿真只能解释一些均匀磨损、不对中的特征真正的现场问题往往出在局部缺陷——外圈滚道剥落、内圈滚道裂纹、滚动体凹坑。建模时有几种流派位移激励法、刚度突变法、附加力脉冲法。初学者适合先从位移激励法入手。位移激励法的物理图像是这样的滚动体滚过缺陷坑时内圈质心会额外下沉一个微小位移相当于在该滚动体的接触路径上叠加一个负向位移偏移。缺陷产生的“撞击”源自几何突变的位移补偿而非额外脉冲力。假设外圈滚道在角度 ( \theta_f ) 处有一个矩形微凹坑深度为 ( d_f )。当第 ( j ) 个滚动体位于缺陷角度附近时它的等效变形量要减去缺陷深度[ \delta_j x \cos \theta_j y \sin \theta_j - c_0 - d_f \cdot I(\theta_j, \theta_f, \Delta\theta) ]其中 ( I ) 是指示函数当角距离落在缺陷宽度 ( \Delta\theta ) 内时取1。实际实现时要注意角度是周期性的直接用 ([\theta_f-\Delta\theta/2, \theta_f\Delta\theta/2]) 判断会出错跨越0度和360度边界时需要先做一个角度归一化差值处理。内圈缺陷处理起来要稍加小心因为内圈本身在旋转缺陷位置随时间以 ( \omega_{shaft} ) 旋转所以缺陷角位置要写成 ( \theta_f \omega_{shaft} t )。而滚动体公转频率是 ( \omega_c )不同频率的耦合会导致内圈故障特征频率 ( BPFI ) 的调制边带——这正是频谱上内圈故障旁边总能看到以转频为间隔的边带的原因。能把这个边带跑出来说明模型物理规律捕捉到位了。4.2 从时域到包络谱验证特征频率故障仿真做完直接在位移信号上看缺陷冲击很不直观。工程上标准的做法是先做包络分析。对加速度信号做窄带滤波后求包络再对包络做FFT得到的频谱在特征频率处会出现明显峰值。在你的仿真代码里可以这样走一套流程对位移信号 ( y(t) ) 二次差分得到加速度然后带通滤波比如500Hz到5000Hz再用Hilbert变换取包络最后FFT。Matlab里这些都有现成函数acc diff(Yq, 2) * fs^2; acc acc - mean(acc); % 设计一个带通滤波器比如二阶巴特沃斯 [b, a] butter(2, [500 5000]/(fs/2), bandpass); acc_f filtfilt(b, a, acc); env abs(hilbert(acc_f)); % 对包络做FFT N length(env); f_axis (0:N-1)/N * fs; Env_fft abs(fft(env));注意diff会把序列缩短2个点做FFT前要对齐长度或者干脆用detrend后的粗估计。包络谱上应该在BPFO处出现明显谱峰。如果实际操作发现特征频率幅值淹没在噪声中一个有效的办法是增加平均次数、延长仿真时间让频率分辨率达到0.5Hz以下对1Hz以内的间隔分辨有决定性作用。下表是我在做某外圈缺陷仿真时的频率对应关系参考特征理论值(Hz)仿真谱峰(Hz)误差说明转频 ( f_r )3030一致外圈故障频率 BPFO约487487峰值对应保持架频率 FTF约11.8511.8低分辨率下易漏滚动体故障频率 BSF约179178.5因滑动假定稍有偏差初学阶段如果看到峰值位置差了一两个频率分辨率格不用太紧张先检查频率分辨率 ( \Delta f 1/T_{sim} )确保仿真时间足够长再做精细比对。5. 常见问题与排查技巧实录5.1 求解器报错和NaN之谜ode45在模型刚度过高、初始瞬态过猛时经常在初始阶段报“步长太小”或产生NaN。原因往往是初始位移设得太大导致接触力瞬时激增。我当时的调试思路是先把初始位移全部置零把重力加载过程人为拆成斜坡前0.02秒重力从0渐变到全值给系统一个“软启动”。这个方法屡试不爽尤其适合新手排查数值刚性问题。如果全程NaN检查接触力函数中是否有除零、指数运算边界溢出等情况。( \delta^{3/2} ) 在 ( \delta ) 为负时如果忘记截断Matlab会得到复数求解器直接崩。可以先在函数里加一句delta max(delta, 0);而不是用 if 包裹整个力计算。5.2 特征频率对不上号仿真出了漂亮的频谱但峰不在理论频率上是新手最容易困惑的事。这时先回头检查保持架转速 ( \omega_c ) 计算公式里的直径比。深沟球轴承的保持架频率是 ( \frac{f_r}{2}(1 - \frac{d_b}{D}) )如果你写的 ( \omega_c ) 是从角速度形式直接代入轴转频就会差一个 ( (1 - d_b/D) ) 的因子。另一个排查点是滚动体个数。理论上如果滚动体均匀分布系统会以 ( n \cdot FTF ) 为基频但实际接触变形非线性会使基频出现在 ( FTF ) 及其整数倍谱图上可能既有细小的 ( FTF ) 峰也有大能量的 ( BPFO ) 峰。这不是代码错了是物理现象。初学者不必为多出来的小峰慌张。5.3 结果总在漂移系统的静态平衡点有时你会看到内圈质心位移一直在往下漂或者整个波形慢慢“沉”下去。这大概率是初始位移离真实平衡点太远且阻尼太小瞬态自由振荡尚未衰减完。增加阻尼系数或延长仿真时间都能改善。这里提醒一句外载荷 ( F_{load} ) 如果远大于重力初始下沉量要用总载荷计算否则平衡点偏差更大。我把几个典型问题整理成速查表方便实际调试时对照现象可能原因解决方案ode45算不动/报错刚度太大或初始瞬态过猛软启动加载或换ode15s结果为NaN变形量未截断导致复数用 max(delta, 0) 截负频谱峰位置偏移保持架转速公式错误核验直径比、转频代入项波形收敛太慢阻尼太小/平衡点偏移增加阻尼比设置初始估计特征频率处没有峰缺陷未正确进入承载区检查缺陷角度与滚动体轨迹5.4 从仿真到实验验证的个人体会我之前用这套模型做外圈缺陷仿真与现场实际采集的信号做对比发现频谱特征频率位置非常吻合但幅值谱的形状有偏差。原因在于实际轴承座传递路径含有多个共振峰而仿真模型中轴承座是刚体不做动力学放大。后期我在内圈质量块和“虚拟传感器”之间加了一组弹簧-质量系统模拟轴承座共振幅值形状就能对上了。个人经验是对于初学者把仿真结果当成“理论基准”比当成“真实数据”更有价值。你可以在仿真上验证信号处理算法包络解调、倒谱分析、深度学习的特征提取确定它能工作后再迁移到现场数据上。省掉了大量盲调算法参数的时间。再给一个执行层面的建议参数扫描时不要一次性跑大批量工况先在单点工况把无缺陷模型调到正常响应再逐步加缺陷深度、改变转速。每一步都保存结果图这样出了问题还能对照是哪一步引入的异常。这套“由简到繁、逐步验证”的思路比我当年直接上完整模型莽撞调试高效得多。最后分享一个小技巧仿真结果如果想导出给报告或后续复现使用最好连同参数结构和版次号一起存成.mat文件。这个习惯帮你省下很多“想不起来当时用什么参数跑出的结果”的尴尬时刻。