干设备振动监测这一行，打交道最多的就是滚动轴承，尤其是深沟球轴承。早期我在某实验室做故障诊断课题时，被一个很现实的问题卡住了：手上有一堆现场采集的轴承振动数据，想验证自己设计的诊断算法，但正常理论解析又没有解析解，拿实验数据一对照又总被各种工况噪声干扰。后来一个偶然机会接触到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微米。初始位移设置为 ( x=0, 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以内的间隔分辨有决定性作用。

下表是我在做某外圈缺陷仿真时的频率对应关系参考：

初学阶段如果看到峰值位置差了一两个频率分辨率格，不用太紧张，先检查频率分辨率 ( \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} ) 如果远大于重力，初始下沉量要用总载荷计算，否则平衡点偏差更大。

我把几个典型问题整理成速查表，方便实际调试时对照：

5.4 从仿真到实验验证的个人体会

我之前用这套模型做外圈缺陷仿真，与现场实际采集的信号做对比，发现频谱特征频率位置非常吻合，但幅值谱的形状有偏差。原因在于实际轴承座传递路径含有多个共振峰，而仿真模型中轴承座是刚体不做动力学放大。后期我在内圈质量块和“虚拟传感器”之间加了一组弹簧-质量系统，模拟轴承座共振，幅值形状就能对上了。

个人经验是，对于初学者，把仿真结果当成“理论基准”比当成“真实数据”更有价值。你可以在仿真上验证信号处理算法（包络解调、倒谱分析、深度学习的特征提取），确定它能工作后，再迁移到现场数据上。省掉了大量盲调算法参数的时间。

再给一个执行层面的建议：参数扫描时不要一次性跑大批量工况，先在单点工况把无缺陷模型调到正常响应，再逐步加缺陷深度、改变转速。每一步都保存结果图，这样出了问题还能对照是哪一步引入的异常。这套“由简到繁、逐步验证”的思路，比我当年直接上完整模型莽撞调试高效得多。

最后分享一个小技巧：仿真结果如果想导出给报告或后续复现使用，最好连同参数结构和版次号一起存成.mat文件。这个习惯帮你省下很多“想不起来当时用什么参数跑出的结果”的尴尬时刻。