简介齿轮动力学关注齿轮传动系统的振动、冲击与噪声机理是机械故障诊断与优化设计的基础。这份MATLAB源码包主要面向机械类专业学生、齿轮传动设计人员及振动信号处理入门者帮助其快速搭建齿轮动力学仿真与短时傅里叶变换STFT分析流程。资源共3个文件包含2个.m脚本和1个.asv自动保存文件整体仅5KB脚本精简以基础模型为主便于逐行理解微分方程求解、齿轮啮合周期激励构造以及STFT加窗取谱的处理步骤可直接在MATLAB中运行修改。资源已有3139人次浏览学习验证了其在教学和科研场景中的参考价值。掌握该源码后读者可进一步分析不同齿数、模数、压力角和负载下的动态响应识别齿面磨损、裂纹、不对中等故障引发的边频带成分并扩展小波变换、多尺度分析等高级方法为设备状态监测与可靠性设计提供依据。1. 齿轮动力学为什么正常齿轮的振动信号里全是“边频带”一个刚出厂、没有任何缺陷的齿轮箱测出来的振动频谱也不是一条干净的啮合频率谱线而是一簇以啮合频率为中心、向两侧展开的“边频带”。很多刚接触齿轮动力学的人第一次看到这张图都以为设备坏了其实这恰恰是齿轮传动的本性——啮合刚度随时间变化、负载波动、转速轻微脉动都会对振动信号产生调幅调相让能量从中心频率向两边扩散。读不懂这个现象后面做故障诊断时就会把正常信号误判成缺陷把真实缺陷当噪声滤掉。这就要引入短时傅里叶变换STFT这个工具。齿轮振动是非平稳信号转速从 800 rpm 拉到 3000 rpm 的过程中啮合频率本身在漂移普通傅里叶变换把整段时间压成一张静态频谱等于把不同时刻的频率信息混在一起平均。STFT 把信号切成短段、逐段做傅里叶变换才能看清频率随时间怎么走。这份 MATLAB 齿轮动力学源码把“动力学模型建模 振动信号生成 STFT 时频分析”串成一条完整链路搞懂它你就同时掌握了齿轮振动的产生机制和它的分析工具。2. 动力学模型与啮合激励刚度波动如何制造振动2.1 从几何参数到“时变啮合刚度”齿轮啮合过程不是一个刚体接触而是一个弹性接触过程。齿面在载荷下发生接触变形轮齿像一个弯曲的悬臂梁同时还有 Hertz 接触变形和齿根圆角处的应力集中。把这些弹性效应折算到啮合线上就得到一个等效弹簧它的刚度不是常数而是随啮合位置周期变化。为什么刚度会波动因为实际齿轮的重合度通常在 1 到 2 之间意味着啮合过程中有时一对齿承担全部载荷有时两对齿同时承担。单齿啮合区刚度低双齿啮合区刚度高这个高低交替就形成了周期性的刚度激励。齿轮旋转一圈每个轮齿经历一次从进入到退出啮合的过程所以刚度波动的基频就是啮合频率fm (z * n) / 60其中 z 是齿数n 是转速单位 rpm。比如一个 23 齿的主动轮在 1500 rpm 下转动啮合频率就是 575 Hz。但这个 575 Hz 不是一条干净的谱线——刚度波形不是标准正弦它接近一个带平顶的周期方波傅里叶展开后自然包含 2 倍、3 倍啮合频率的高次谐波。这就是为什么齿轮振动频谱里经常能看到二倍频、三倍频的原因不必把它们一律当成故障特征。2.2 单自由度扭转模型把齿轮副抽象成什么理解齿轮动力学不需要一上来就上有限元。工程上最常用的初步模型是单自由度扭转振动系统把主动轮和从动轮看成两个转动惯量中间通过时变啮合刚度相连再并联一个阻尼器。传动误差包括齿面形貌误差、装配误差、轮齿变形作为位移激励加到系统的右侧。J1 * θ1 c * (θ1 - θ2) k(t) * (θ1 - θ2) T1 J2 * θ2 - c * (θ1 - θ2) - k(t) * (θ1 - θ2) -T2把两个方程相减定义一个相对位移变量就能化成一个等效的单自由度方程me * x c * x k(t) * x F(t)这里 me 是等效质量k(t) 是时变啮合刚度F(t) 是综合激励。用 MATLAB 的 ode45 可以直接数值求解这个变系数微分方程。这个模型虽然简单却抓住了齿轮振动最核心的东西刚度激励加上误差激励。齿面磨损、点蚀、裂纹最终都表现为刚度波形或误差波形的改变因此这个模型是后续一切故障诊断的出发点。2.3 边频带的形成机制调幅与调相的合力回到开头那个现象为什么正常齿轮会有一簇边频带答案在于啮合刚度不仅按啮合频率变化还同时被轴频调制。转轴每转一圈由于偏心、质量不平衡、齿距偏差的存在啮合力的幅值会轻微起伏。这个起伏的频率是轴频 fa n / 60它像一个低频包络骑在啮频载波上这就是标准的调幅模型y(t) (1 m * cos(2π * fa * t)) * cos(2π * fm * t)对调幅信号做傅里叶变换频域里会出现三根谱线中心频率 fm两侧各一根 fa 间隔的边带。齿轮的刚度波动里还包含相位调制成分于是边带进一步丰富形成以 fm 为中心、以 fa 为间隔的边带族。这就是齿轮振动频谱“一簇一簇”的由来。理解了这一点故障诊断的基本逻辑就清楚了边带间隔等于哪个轴的转频故障就在哪个轴上边带幅值相对于中心频率的比值变大说明调制深度变大往往意味着齿面缺陷在发展。只看幅值不看边带间隔的频谱分析在齿轮故障诊断里基本是盲人摸象。3. 傅里叶到短时傅里叶为什么常规频谱图看不清故障3.1 傅里叶变换的全局平均缺陷普通傅里叶变换把整个时域信号分解成不同频率的正弦波之和得到的是一个“时间无关”的频谱。这个性质决定了它只适合平稳信号——频率成分不随时间变化的信号。齿轮箱升降速过程中转速连续变化啮合频率从 500 Hz 一路升到 1500 Hz变频过程中还可能短暂共振。对整段信号做一次 FFT结果是把不同时刻的频率成分叠加到同一张图上任何瞬态冲击、共振穿越、转速漂移都被平均掉了。更麻烦的是齿轮故障信号本身就是瞬态的。齿面裂纹在啮合的某个瞬间产生一次冲击这个冲击持续不到几毫秒能量分散在整个频带里。全局 FFT 会把这次冲击的能量摊薄到整段时间上导致它在频谱里几乎看不见。这就是为什么有经验的诊断工程师从来不看整段信号的静态频谱而是看时频图或转速跟踪谱。3.2 加窗切割短时傅里叶变换如何工作短时傅里叶变换的思路很直接假设信号在足够短的时间窗内近似平稳用窗函数从信号里截出一小段对这一段做 FFT得到这一时刻附近的频谱。然后滑动窗函数重复操作把每一帧的频谱按时间排列就得到一张“频率-时间-幅值”三维谱图也就是瀑布图或时频图。S(t, f) ∫ x(τ) * w(τ - t) * e^(-j2πfτ) dτ窗函数 w(τ - t) 把关注点收窄到 t 附近。MATLAB 的 spectrogram 函数把这一过程封装成一步调用用户在大多数情况下只需要指定窗长、重叠比例和 FFT 点数。窗长直接决定时间分辨率和频率分辨率的取舍窗越长频率分辨率越高能分辨靠得很近的边频带但时间窗变长频率随时间的变化会被窗口内的平均模糊掉时间分辨率变差。这两个分辨率不可能同时提高因为海森堡不确定性原理在信号处理里的投影就是时间带宽积存在下限。3.3 窗函数家族Hamming、Hanning、Blackman 到底改了什么选择窗函数本质上是选择频谱泄漏的形态。矩形窗旁瓣衰减只有 -13 dB主瓣窄但旁瓣高强频率成分的旁瓣会掩盖旁边的弱信号。Hamming 窗把旁瓣压低到约 -43 dBHanning 窗略低一点而 Blackman 窗可以压到 -58 dB。代价是主瓣变宽本来能分开的边带可能糊在一起。齿轮振动分析中常见的做法是先用 Hanning 窗观察整体时频结构因为它旁瓣低、主瓣宽度适中适合看变频过程中的啮频轨迹如果要做精细的边带间隔测量再换更长的矩形或 Hamming 窗提高频域锐度。同一份数据用不同窗长看出来的结论可能完全不一样——窗长选太短啮频轨迹在时频图上是粗壮的条带边带看不清窗长选太长升降速过程中相邻两秒的频率变化被平均掉时频图变得模糊。4. MATLAB 源码剖析chilundonglixue.m 里每一步在算什么4.1 文件结构与自动化备份的启示这份资源的核心文件是 chilundonglixue.m 和 chilundonglixue.asv。.asv 是 MATLAB 编辑器在自动保存时生成的备份文件说明代码经历过较长周期的迭代调试——这不是一次性写完就发布的脚本而是边跑边改的产物。还有一个以 428ec93e6a72e45941d3c19a586373e7.m 命名的文件这种长哈希文件名通常是 MATLAB Online 或版本控制中的临时文件把它和主脚本对照着看能看到作者在模型参数调整过程中的留痕。主脚本的功能链路分四段定义齿轮和工况参数、构造时变啮合刚度、用数值积分求响应、对响应做 STFT 可视化。每一段对应齿轮动力学分析的一个标准步骤。4.2 参数定义与激励构建% 齿轮基本参数 z1 23; % 主动轮齿数 z2 47; % 从动轮齿数 n1 1500; % 主动轮转速, rpm fa n1 / 60; % 轴频 fm (z1 * n1) / 60; % 啮合频率 % 时变啮合刚度近似单双齿交替形成方波正弦修正 k0 2e8; % 平均啮合刚度, N/m epsilon 1.5; % 重合度 k k0 * (1 0.3 * sin(2*pi*fm*t pi/2));这里啮合刚度用一个在平均值附近波动的正弦近似虽然幅度和相位做了简化但抓住了刚度激励“以啮合频率周期变化”的本质。重合度对刚度波动幅度有直接影响重合度越接近整数刚度波动越大振动越剧烈重合度越远离整数双齿啮合占比越高刚度越平稳。4.3 ode45 求解与振动响应% 等效单自由度振动系统 % x 2*zeta*omega_n*x omega_n^2 * x F(t)/me me 0.5; % 等效质量估算 zeta 0.03; % 阻尼比 omega_n sqrt(k0 / me); % 固有频率 Fn 0.5 * k0 * (1 - cos(2*pi*fm*t)); % 啮合冲击力 [tau, Y] ode45((t, y) gear_eq(t, y, fm, me, zeta, omega_n), ... [0 1.0], [0 0]);ode45 是变步长积分器误差控制由相对容差和绝对容差决定。齿轮方程里刚度以高频振荡形式存在步长必须足够小才能捕捉到啮合频率的完整周期。如果发现响应波形在啮合频率附近出现锯齿状跳动第一反应不是改阻尼参数而是把求解器的最大步长 MaxStep 调小到啮合周期的十分之一。4.4 spectrogram 时频分析的实际调用% 短时傅里叶变换 fs 10000; % 采样率 10 kHz window hamming(256); % 每帧 256 点约 25.6 ms noverlap 220; % 85% 重叠提升时间分辨率 nfft 1024; % FFT 补零到 1024 点细化频域网格 [S, f, t_stft] spectrogram(x, window, noverlap, nfft, fs); imagesc(t_stft, f, 20*log10(abs(S)eps));noverlap 取 220 表示相邻窗只错开 36 个采样点时间轴步长由 Hop Size 决定与 FFT 点数无关。补零到 1024 点不会提高真实物理分辨率只是让频谱看起来更“顺滑”。真正决定是否能分辨边带的是窗长256 点在 10 kHz 采样率下对应频率分辨率约 39 Hz如果轴频只有 25 Hz边带间隔小于频率分辨率边带就融进主瓣里看不见了。这时要么加长窗到 512 或 1024 点要么降采样提高每单位时间内的有效点数。5. 避坑清单齿轮振动分析里常见的四个翻车现场5.1 窗长缩短之后边带神秘消失现象同一组数据把窗长从 1024 改成 256 后频谱上的边频带突然不见了画面干净得反常。有人拿着这张图去报告“齿轮状态良好”。原因窗长 256 点对应的频率分辨率只有约 39 Hz而边带间隔等于轴频 25 Hz两个边带和中心频率的间距都小于分辨率三根谱线被主瓣合并成一根视觉上边带被“抹平”了。解决先算清轴频和期望频率分辨率的匹配关系。频率分辨率约等于采样率除以窗长要让边带可见窗长至少要满足“分辨率小于 0.5 倍边带间隔”。诊断齿轮故障时优先保证频域分辨率时间分辨率居其次。5.2 变速工况下时频图出现“假共振”现象设备从 800 rpm 匀加速到 3000 rpm时频图上出现一条明亮的弧线看起来像是某阶固有频率被激发了。原因这条弧线大概率不是共振而是啮合频率的斜坡轨迹。匀加速过程中啮合频率随时间线性上升在时频图上天然就是一条斜线。它的亮度与 FFT 幅值成正比容易被误读为共振带。解决把理论啮合频率轨迹叠加到时频图上做对照如果亮带与理论轨迹重合而固有频率附近没有独立亮斑那就是正常扫频而非共振。更严谨的做法是提取瞬时转速信号做阶次跟踪把时频图转成阶次谱再判断。5.3 滚动轴承故障被齿轮边带掩盖现象齿轮箱里同时有轴承损伤和齿轮磨损频谱图上齿轮的边带族幅值明显轴承故障特征频率完全被淹没。原因轴承故障信号的幅值通常比齿轮啮合信号低一个数量级而且传递路径远、衰耗大普通频谱里能量集中在低频段。解决对信号做包络分析先带通滤波到轴承特征频率所在频段再做包络谱或 STFT。包络解调把高频冲击解调成低频包络相当于把微弱故障信号从“载波”里提取出来这是轴承诊断的标准操作配合时频图可以区分齿轮故障在高频段还是轴承故障的高频冲击特征。5.4 相位信息被丢弃导致时域盲区现象只看幅频谱和时频图判断故障严重程度时发现幅值变化不大但设备噪声明显加剧。原因FFT 和 spectrogram 输出的是幅值谱或功率谱相位信息全部丢弃。齿轮早期裂纹有时首先改变的是相位调制特征即频率相同但相位随时间偏移幅值频谱上反应微弱。解决把原始时域信号分段做重采样对齐计算各段之间的相位差。某一段相位出现跳变或漂移往往是局部缺陷的早期征兆。保持在分析流程里同时输出时域包络和相位曲线不要只盯幅值图。6. 设计取舍用时间分辨率换频率分辨率这笔交易怎么才划算STFT 的窗长选择本质上是拿时间精度换频率精度这是它的数学模型决定的没有免费午餐。海森堡不确定性原理在这里表现为时间窗标准差和频率窗标准差的乘积不小于某个常数你不可能让两者同时变小。问题是齿轮诊断场景里如何判定交易划不划算。判断标准取决于你要抓的目标。做升降速过程监测关心的是啮频轨迹的跟随性这时候时间精度优先窗长取短哪怕频率分辨率差一点只要能看到斜线连续不断就行。做定转速下的边带分析关心的是每 25 Hz 一根边带是否清晰分离这时频率精度优先窗长加长牺牲时间维度上的瞬时性。最怕的是不分场景固定一套参数拿到什么数据都用同一窗长跑结果两头不讨好。我一般会这么定流程先用短窗跑一次全局时频图确认信号里有哪些频率成分、大致的时间分布再用长窗对感兴趣的时间段做局部精细分析测边带间隔。两次分析用同一份数据窗长可能差 4 到 8 倍但得到的信息互补。这套“全局粗看、局部细看”的打法比纠结某一个窗长参数要靠谱得多。如果想在这个资源基础上进一步延伸下一步值得做的是替换小波变换。STFT 的窗长固定意味着高频段和低频段的时间分辨率一样但齿轮故障冲击是宽频信号——低频段需要长窗保频率精度高频段需要短窗定位冲击时刻这是一个矛盾。小波变换用尺度参数替代窗长在高频处自动用短窗、低频处自动用长窗恰好缓解这个矛盾。用 MATLAB 的 cwt 函数对比同一组振动数据你能直观看出 STFT 在某段频率上的模糊现象被小波改善到什么程度。从那以后我每次处理齿轮数据都强制先跑一遍 STFT 全局图再决定局部参数已经成了习惯希望帮到你。本文还有配套的精品资源点击获取