MATLAB滚动轴承动力学仿真详解：正常/外圈/内圈/滚动体故障建模与验证

做滚动轴承故障诊断的人，手头最缺的东西往往不是数据，而是一套能讲清楚“故障信号到底怎么来的”仿真模型。很多刚接触这个方向的朋友上来就找公开数据集，下载一堆实测信号开始跑包络谱，跑出来对不上故障特征频率，也不知道是数据问题还是算法问题。我在这个方向摸爬滚打了好几年，最大的体会就是：先有一套可靠的仿真代码，你才知道真实信号里哪些特征是本质的、哪些是噪声带来的。这篇博文要分享的，就是一套基于MATLAB的滚动轴承动力学仿真代码，覆盖正常、外圈故障、内圈故障、滚动体故障四种状态，完全按照滚动轴承故障机理建模，不是简单地用MATLAB生成几个正弦波叠加就完事。代码的定位很明确：给做故障诊断、状态监测、信号处理算法验证的工程师和学生，提供一个“已知答案”的信号源——每类故障的故障特征频率都能从机理上解析计算出来，仿真结果可以直接对照理论值验证，也能用来测试包络谱、时频分析、深度学习分类等后续算法的效果。

1. 轴承故障模拟的整体思路与设计拆解

1.1 为什么一定要从故障机理出发搭建仿真模型

先聊一个很现实的问题：很多人觉得轴承故障仿真就是按照故障特征频率公式，生成几个正弦衰减振荡叠加在一起，再加点白噪声就完事了。这个思路确实能跑，但做出来的信号太“干净”了，跟实际采集的振动信号脱节严重。实测信号里那种非平稳、调幅调频、随机滑移的复杂表现，是简单叠加模型根本复现不出来的。

我做的这套代码，出发点是把轴承看作一个机械动力学系统：内圈固定在转轴上随轴一起旋转，外圈安装在轴承座或箱体上保持静止，滚动体在内外滚道之间既公转又自转。当某个元件出现局部损伤时，滚动体滚过损伤位置会产生一个周期性的冲击力，这个冲击力会激起轴承系统自身的固有振动。不同的故障位置，冲击力的作用方式和路径不同，最终在传感器位置拾取到的振动信号自然也就有本质差异。把这些差异从数学和力学层面建模出来，才是“根据故障机理建模”的真正含义。

从实际工程角度看，这套基于机理的仿真有一个巨大优势：可解释性。每个参数都有明确的物理意义——轴承几何尺寸、转速、载荷、损伤尺寸和位置。你用代码跑出来的信号，理论故障特征频率是可以通过公式精确计算的，后面做包络分析、机器学习分类，都有“标准答案”可以对照。

1.2 支撑这套模型的三个基础理论

要真正看明白代码，有三个理论基础必须理清，这三个基础理论也是整个建模方案的骨架。

第一个是滚动轴承的运动学关系。滚动体在内外滚道之间运动时，接触点上的速度关系决定了滚动体公转速度，也就是保持架的转速。这个转速和轴承各元件的故障特征频率直接挂钩——外圈固定时，外圈故障特征频率等于保持架转速乘以滚动体个数；内圈故障时，因为内圈在旋转，特征频率要再加一个转频项。这些关系是经典文献里的标准公式，代码里全部用几何参数直接计算。

第二个是冲击响应模型。每个故障产生的冲击不是理想的脉冲，而是激发出轴承系统的高频固有振荡，这个振荡会按指数规律衰减。用一句话概括：故障冲击是激励，系统固有频率是载体，传感器拾取到的是两者的卷积结果。因此代码里需要设定系统的固有频率和阻尼比，仿真出来的信号才是“有血有肉”的衰减振荡，而不是孤零零的尖峰。

第三个是载荷分布与传递路径的影响。轴承在工作时承受径向载荷，只有位于载荷区内的滚动体才能产生有效的接触力。外圈固定时，故障点位置固定，但载荷作用下每个滚动体进入载荷区的受力大小不一样；内圈旋转时，故障点跟着内圈转，形成了“定期通过载荷区”的调幅现象。这一层影响决定了信号的幅值包络特征，是区分内圈/外圈故障的关键物理因素，也是这套代码里比较注重的建模细节。

1.3 模型选型：为什么选集中参数模型而不是有限元

一开始我也纠结过：要不要用ANSYS或者Abaqus做精细的有限元模型，把轴承每一个接触应力都算出来？后来实践证明，对故障诊断信号仿真而言，有限元的成本完全不成比例。一套足够精细的轴承-转子有限元模型，建模周期以周为单位，单次瞬态仿真可能要跑几个小时，而且精度受接触算法、网格质量影响非常大。

相比之下，集中参数模型（lumped parameter model）用等效质量、刚度、阻尼来描述轴承系统的动态特性，计算量极小——MATLAB里跑一次几秒钟就能完成几千转的数据。它牺牲的是微观应力分布的精度，但完整保留了故障冲击在振动信号中的核心特征：周期性、衰减振荡、调制效应。对故障诊断研究来说，我们要的本质就是这些特征，而不是某个接触点上的应力数值。

代码中采用的模型本质上是把轴承座和传感器拾取路径等效成一个或几个二阶振荡系统，故障冲击作为激励输入。这个方案在学术论文和工程实践中被广泛验证过，不仅能快速生成足够真实的仿真信号，还能灵活调节各种参数来模拟不同工况，这才是故障诊断算法研究需要的“信号发生器”。

2. 四种轴承状态的建模差异与核心细节

2.1 正常状态：不是简单的一条直线

健康的滚动轴承在运转时，即使没有任何局部缺陷，振动信号也不是零。我在代码里把正常状态分成三个分量来建模，这比很多人想的要复杂一层。

第一是转频及其谐波分量。转子系统不可避免地存在残余不平衡，这会产生以转频为基频的振动，还会伴随少量二次、三次谐波。这个分量在频谱上非常容易识别，频率等于转速除以60，单位是Hz。第二是滚动体通过频率分量。即使滚道表面光滑，滚动体通过时也会由于有限的滚动体个数而产生周期性的刚度变化，这个分量的频率跟滚动体个数和保持架转速有关，幅值通常比较小。第三是随机噪声和低频漂移。轴承座、基础和其他机械部件的干扰，以及实测环境中的随机波动，用白噪声加低频振荡来模拟。

正常状态仿真最重要的价值在于提供一个“基准线”。后面做故障诊断，所有的特征提取和阈值判断都要拿正常状态对照。代码里正常状态的输出信号，会保持幅值平稳、频谱成分单纯，这是后续三种故障状态共同的“底噪”。

2.2 外圈故障：固定损伤点的冲击序列

外圈故障（BPFO，Ball Pass Frequency Outer Race）的机理在四种状态里最直接。外圈安装在轴承座上，位置固定，所以故障点相对于传感器的位置基本不变。滚动体每经过一次故障点，就产生一次冲击，这个冲击的间隔时间完全相等，对应频率就是外圈故障特征频率：

• BPFO = (n / 2) × fr × (1 - (d / D) × cos α)

其中n是滚动体个数，fr是轴转频，d是滚动体直径，D是节圆直径，α是接触角。

代码在建模外圈故障时，主要注意两点：一是故障点位置落在载荷区之内，因为只有载荷区的滚动体才能形成足够的接触力，所以冲击幅值基本平稳；二是冲击间隔非常均匀，因为外圈不转动，一旦滚动体转过一周，经过故障点的时间间隔就是严格周期性的。仿真出的信号在时域上是等间隔的衰减冲击串，包络谱上会在BPFO处有清晰的谱峰，还会伴随BPFO的二倍频、三倍频谐波。这个规律性，正是实测外圈故障中比较容易识别、也最容易验证的。

2.3 内圈故障：旋转载荷区的调制效应

内圈故障（BPFI，Ball Pass Frequency Inner Race）的机理要比外圈复杂不少。内圈随轴一起旋转，故障点也在旋转，但径向载荷方向是固定的，轴承内部的载荷分布也随位置固定。这样一来，故障点每次转动到载荷区时，冲击比较强；转到非载荷区时，冲击明显变弱甚至接近消失。这个强弱交替的过程，在信号上就体现为“载荷调制”。

因此内圈故障的冲击序列属于典型的调幅信号，载波是衰减振荡的高频固有振动，调制信号是转频相关的周期性包络变化。在频谱上，除了BPFI这个核心频率之外，BPFI两侧会出现以转频fr为间隔的边带族，形成“f_BPFI ± k×fr”的谱线结构。

代码对内圈故障的处理，是在每次冲击的幅值上乘以一个与故障点位置相关的权函数。这个权函数不是简单的正弦波，而是根据载荷分布形状（通常是中间强两端弱的类余弦形状）来构造。这样仿真出来的信号，包络谱上不但能看见BPFI，边带的分布也跟实际测得的内圈故障信号特征对应得比较准。

2.4 滚动体故障：动圈耦合与双冲击

滚动体故障（BSF，Ball Spin Frequency）在三类故障里最容易搞混，因为它牵涉到滚动体的自转。滚动体自身的缺陷随滚动体一边公转一边自转，故障点并不始终与滚道接触。

• BSF = (D / (2d)) × fr × (1 - (d / D)² × cos² α)

滚动体故障的名义故障特征频率是BSF，但实际信号里有两个格外值得关注的现象。

第一个现象是故障位置的周期性暴露：滚动体自转一周，故障点会各接触一次内滚道和外滚道，所以在理论上每次旋转会产生两次冲击，实际冲击的周期其实跟滚动体自转周期一致，而BSF是滚动体自转频率的2倍关系。所以代码里的冲击频率要用自转频率而非BSF。

第二个现象是冲击幅值的缓变调制：滚动体在公转过程中会进出载荷区，每次进入载荷区时，与外滚道和与内滚道的接触力都比较大，离开载荷区后冲击显著减弱。这就是为什么滚动体故障信号里还叠加了一层以保持架转速为周期的慢调制。代码里这两层关系都要体现，仿真信号才能找到那个既不是BSF也不是转频、而是二者组合的复杂频率结构。这也是很多初学者觉得滚动体故障仿真最难的地方。

3. MATLAB代码实现与关键参数设置

3.1 整体程序结构与函数划分

整套代码我按模块化思路来写，每个模块独立成一个函数文件，方便单独调试和替换。整个工程的主程序结构如下：

%% 主程序：轴承动力学仿真入口 % 步骤1：定义轴承几何参数和工况参数 % 步骤2：调用函数计算故障特征频率 % 步骤3：根据状态类型生成冲击序列 % 步骤4：调用冲击响应函数卷积生成振动信号 % 步骤5：添加噪声、趋势项，输出时域波形 % 步骤6：调用分析函数做包络谱或频谱验证

代码文件划分如表所示：

第一次上手的人，我建议直接跑主函数simulate_bearing_signal.m，先看输出结果，再一头扎进子函数里理解细节。别一开始就试图把每个模块都弄透，先建立整体感知，再深挖原理。

3.2 滚动轴承几何参数与故障特征频率计算代码

故障特征频率是整个仿真信号的灵魂，代码里用到的都是滚动轴承运动学里最基础的计算公式。我用一个内置的深沟球轴承参数作为默认示例，以大家常用的6205轴承为例：

function [BPFO, BPFI, BSF, FTF] = calc_fault_freqs(d, D, n, alpha, fr) % 输入： % d - 滚动体直径 (mm) % D - 节圆直径 (mm) % n - 滚动体个数 % alpha - 接触角 (rad) % fr - 轴转频 (Hz) % 输出： % BPFO - 外圈故障特征频率 (Hz) % BPFI - 内圈故障特征频率 (Hz) % BSF - 滚动体故障特征频率 (Hz) % FTF - 保持架故障特征频率 (Hz) c = cos(alpha); BPFO = n / 2 * fr * (1 - d / D * c); % 外圈固定时，滚动体通过外圈频率 BPFI = n / 2 * fr * (1 + d / D * c); % 内圈旋转时，滚动体通过内圈频率 BSF = D / (2 * d) * fr * (1 - (d / D * c)^2); % 滚动体自转频率 FTF = 1 / 2 * fr * (1 - d / D * c); % 保持架公转频率 end

用6205轴承的典型参数（d=7.94mm，D=39.04mm，n=9，α=0，fr=29.3Hz，对应约1758rpm）代入，算出来的结果是：BPFO约105.6Hz，BPFI约158.1Hz，BSF约68.9Hz。这三个频率就是后面判断仿真信号“对不对”的标尺。

注意：故障特征频率公式的前提是外圈固定、内圈旋转，运行时滚动体与滚道之间为纯滚动接触。实际中轻微打滑会导致实测频率与理论值有零点几赫兹的偏差，代码里用随机滑移来模拟这一现象。

3.3 冲击序列生成的核心逻辑

这一步是整个代码里最关键也最容易理解错的地方。我的做法是：先生成一个时间轴，然后判断“哪些时刻该有冲击”，在对应时刻放置冲击脉冲。但这里有三个细节必须处理，否则代码就是“玩具”。

第一个细节是随机滑移。真实轴承中滚动体在滚道上不是纯滚动，存在微小的随机滑移，导致冲击间隔不是绝对均匀，而是在理论周期附近有小幅波动。代码里通过给冲击的到达时间加一个随机扰动来模拟，扰动幅值一般取理论周期的0.5%到2%。

第二个细节是冲击不是瞬时完成的。故障冲击本身有一个持续时间，这个时间跟滚动体滚过故障缺口的长度和速度有关。代码里把单个冲击建模成一个极短的高斯脉冲或半正弦脉冲，而不是一个理想的迪拉克脉冲。

第三个细节是故障相位偏移。尤其是内圈故障，初始相位决定了第一个冲击发生在什么位置，从而影响后续所有冲击的包络形态。代码里设置了故障初始角这个输入参数，用户可以通过改变它观察不同相位下仿真信号包络谱的微小变化。

下面给出一个简化的冲击序列生成代码片段：

function impulse_train = gen_impulse_train(t, fs, fault_freq, amp, slip_ratio, phase) % t: 时间轴 % fs: 采样率 % fault_freq: 故障特征频率（外圈/内圈/滚动体对应各自频率） % amp: 名义冲击幅值 % slip_ratio: 随机滑移比例，例如0.01表示1%滑移 % phase: 初始相位 impulse_train = zeros(size(t)); T_fault = 1 / fault_freq; % 理论冲击周期 k = 0; while k * T_fault <= t(end) % 引入随机滑移 jitter = slip_ratio * T_fault * randn(1); t_impulse = k * T_fault + jitter + phase / (2*pi*fault_freq); if t_impulse >= t(1) && t_impulse <= t(end) [~, idx] = min(abs(t - t_impulse)); % 找到最近的时间索引 % 用一个短高斯脉冲模拟单个冲击的时域形状 tau = 0.0001; % 冲击持续时间（秒） impulse_train(idx) = impulse_train(idx) + amp / (sqrt(2*pi)*tau) * exp(-((t(idx)-t_impulse).^2)/(2*tau^2)); end k = k + 1; end end

3.4 系统脉冲响应与载荷调制的融合

冲击序列生成之后，下一步就是把这一串脉冲变成“像振动信号”的东西。这一步我用了一个关键假设：轴承-轴承座系统可以简化为一个单自由度二阶系统。当一个冲击作用在这个系统上时，输出是一个衰减的正弦振荡：

function h = system_response(fn, zeta, t) % fn: 系统固有频率 (Hz) % zeta: 阻尼比 % t: 时间轴 wd = fn * sqrt(1 - zeta^2); % 阻尼固有角频率 h = exp(-zeta * 2*pi*fn * t) .* sin(2*pi*wd * t); % 将t < 0的部分清零 h(t < 0) = 0; end

系统固有频率fn一般取值在2000Hz到8000Hz之间，因为轴承故障冲击激起的主要是轴承座和传感器的共振频带。阻尼比zeta取值在0.02到0.1之间，阻尼太小，冲击振荡拖得太长，掩盖后续冲击；阻尼太大，振荡衰减太快，信号就变成一串孤立的毛刺。默认推荐fn=3000Hz，zeta=0.05，这个组合跑出来的信号无论时域还是频域都跟实测信号比较接近。

冲击序列与系统响应做卷积，就是整段振动信号的基础。然后对冲击的幅值乘上载荷调制函数——外圈故障用接近常数的权值，内圈故障用载荷区形状的权值，滚动体故障用双周期叠加调制。这部分在代码中体现为：

% 对外圈故障：幅值权值基本恒定，只有轻微波动 % 对内圈故障：权值 = 载荷分布(故障点相对旋转角度) % 对滚动体故障：权值 = 载荷分布 × 内外滚道接触交替系数

做完这些，再叠加正常状态的三个分量，最后按信噪比要求加入白噪声。这样生成的一整段信号，才具备“实测感”。

4. 仿真结果分析与故障特征验证

4.1 时域波形与包络谱的对照判断

跑完仿真之后，第一件事不是直接扔给算法，而是先“肉眼检视”一下信号对不对。我会用两幅图来验证——时域波形和包络谱。这两幅图是判断建模是否合理最快的方式。

时域波形上看三件事：第一，冲击是否周期出现，间隔是否大致稳定；第二，冲击幅值的包络是否符合对应故障类型的调制规律；第三，衰减振荡的形状是否清晰，有没有被噪声完全淹没。如果时域上冲击串都看不清，那后面所有特征提取都不可靠。

包络谱上看两件事：第一，故障特征频率处是否有明显谱峰；第二，谱峰两侧是否出现了理论预期的边带结构。比如内圈故障，若BPFI两侧没有转频间隔的边带，说明载荷调制没加对；外圈故障如果BPFO处只有基频而二倍频三倍频太弱，说明冲击信号的谐波成分不够丰富。

4.2 四种状态的仿真信号特征对比

为了直观对比四种状态差异，我把仿真结果整理成一张特征对照表，这是一次调参过程中记录下来的典型数据（6205轴承，转频29.3Hz，采样率50kHz）：

从表中能看到，每种状态的故障特征频率都跟理论值对得上。内圈故障的边带宽度比外圈故障明显更宽，这是由旋转载荷区的调制作用引起的；滚动体故障的幅值起伏周期明显长于冲击间隔，这对应滚动体进出载荷区的公转周期。如果仿真信号没有体现这些差异，那多半是载荷调制模块的参数没调对。

4.3 仿真信号的可信度验证方法

仿真代码写完，怎么证明它可信？我通常用三个层次的验证方法，建议你也按这个顺序走一遍。

第一层是理论频率对照。直接对比包络谱峰值的最大谱线频率和理论计算值，偏差应该控制在1%以内（这里不包含随机滑移引入的微小偏移）。这一步在代码里写一个自动检测函数，给定仿真信号和理论故障频率，输出峰值附近偏差点。

第二层是周期统计对照。直接从仿真信号的时域波形里，用峰值检测算法提取冲击间隔，再换算成脉冲频率。这个方法可以避开包络谱的分辨率限制，验证信号时长较短、频率较低时依然可靠。

第三层是算法反向验证。把仿真信号输入到独立的包络分析程序或者时频分析工具里，看输出的故障特征频率是否仍能正确识别。这一步实际上是拿仿真信号来验证你自己的分析算法——如果算法连“标准答案”都提不出来，那对实测信号就更不靠谱。

做完这三层验证，这套仿真信号才能说真正可用于后续算法开发了。

5. 常见问题与调试经验总结

5.1 仿真信号看起来不对时的排查顺序

我遇到过不少朋友拿着这套代码来问“为什么我跑出来的信号没有冲击”。排查顺序非常固定，对着这个清单从上往下查，基本能解决90%的问题。

第一，查故障特征频率是否算对了。常见错误是轴承几何参数搞错，比如把滚珠直径填成半径，或者忘记把直径单位换算成米制，导致算出来的BPFO差了整整一倍。

第二，查冲击间隔与总时长的关系。如果仿真时长太短，比如只有0.1秒，外圈故障特征频率105.6Hz，整个信号里只有大约10个冲击，时域上看稀稀拉拉，后面做包络谱分辨率也不够。建议至少仿真1秒以上数据。

第三，查采样率是否足够。采样率必须大于系统固有频率的5到10倍，否则衰减振荡的波形完全失真。代码默认50kHz，实测中一些采集设备是20kHz或10kHz，如果固有频率设在8000Hz，20kHz采样率勉强够用，但10kHz就会明显混叠。

第四，查信噪比参数。噪声太大冲击被淹没，太小信号太干净不像实测。建议从10dB开始试，逐步降低到0dB，看什么信噪比下包络谱还能识别出故障频率。

5.2 与实测信号对不上时的调整方向

做仿真的人最头疼的就是“仿真信号和实测信号长得很不一样”。我的观点是：仿真追求的不是形状上像素级一致，而是统计特征和故障特征的一致性。如果实测信号里BPFO处明显有峰值但旁边还有大量杂峰，仿真信号里BPFO处很干净，这不能说明仿真不对，只能说明实测里的杂峰来自其他干扰源。

如果一定要让仿真更贴近某台具体设备，我建议先调整三个参数：

• 系统固有频率：实测信号中观察共振峰的位置，修改fn到对应频带上，能让信号的整体频域分布更像。
• 阻尼比：阻尼决定实际采集数据中冲击衰减的快慢，从小幅增加阻尼能让信号更“糊”一点，更接近实测的粗糙感。
• 随机滑移比例：实测中接触微滑移导致边带不干净，把滑移比例从0.5%提高到2%，边带会明显展宽，更贴近真实。

5.3 封装函数时的几个实用细节

最后分享几个使用这套代码时的操作细节，这些细节平时写代码不会注意到，但实际跑仿真时能省下不少时间。

第一个细节是固定随机种子。为了复现实验结果，在仿真入口加上rng(2024)这样的语句，保证每次跑同一组参数生成的信号完全一致。这对写论文、对比不同算法非常关键，不然结果永远没法复现。

第二个细节是故障特征频率传参最好不要直接在内部重算。让主函数调用一次calc_fault_freqs算出各频率值，再传给后续的冲击序列生成函数。这样调试时只需要打印一次特征频率，就能确认整条链路用的是同一个频率。

第三个细节是批量仿真建议做成循环加数组存储。把不同故障类型、不同转速的仿真信号存成struct数组或者直接存成.mat文件，后面做深度学习训练时一次性加载，不用临时重跑。

第四个细节是包络谱分析时注意频带选择。仿真信号用的固有频率在3000Hz附近，做包络分析时带通滤波范围建议设在2000到5000Hz之间，这样能最大化保留故障信息，同时避开转频等低频强干扰成分。

我个人在这套代码上调试得最久的地方，其实是滚动体故障的双冲击建模。反复查阅文献、对比实测数据后才发现，单纯用BSF去驱动冲击序列是不对的，必须把滚动体自转频率和公转周期两层因素同时考虑进去。这个坑踩完之后，我对所有故障状态的建模理解都更深了一层。如果你也在这块卡住了，不妨停下来想一想：故障点在哪里，它是怎么运动的，载荷区在哪、变不变——这三个问题想透了，绝大部分建模问题都能解决。