简介面向机械工程与故障诊断领域的MATLAB源码包聚焦齿轮动力学建模与振动信号分析适合需要借助短时傅里叶变换STFT研究齿轮啮合特性、识别齿面磨损或裂纹等故障模式的研究人员和工程师。压缩包共三个文件包括两个.m脚本与一个.asv自动备份文件整体仅5KB体量精巧便于直接阅读和修改。已有3137人学习下载。内容围绕齿轮动力学模型展开通过MATLAB信号处理工具箱中的spectrogram函数实现加窗、切割与可视化可帮助用户掌握从齿轮振动数据中提取时频特征、定位瞬态故障频率的分析思路同时结合窗函数选择、频率分辨率权衡等细节为后续扩展小波变换或有限元分析提供了基础。这套源码对于入门齿轮动力学数值实验和快速搭建振动信号处理流程具有实用参考价值。1. 齿轮动力学建模从一对齿的啮合到整条传动链的振动响应齿轮箱振动超标、减速机异响、齿面早期点蚀这些现场问题背后几乎都指向同一个源头——啮合副的动态激励。手算齿轮强度只能给出静载安全系数回答不了“为什么这台箱体在某个转速下噪声特别大”这类问题。要回答就得建立齿轮动力学模型把齿轮副的质量、刚度、阻尼和激励放到一个微分方程组里用数值方法求解时域响应再从响应里找共振点、边频带和冲击特征。MATLAB在这条路上几乎是默认工具ode45做数值积分、FFT做频谱分析、优化工具箱做参数识别一个环境全链打通。这篇文章按常见做法讲一套能落地的扭转振动模型理论怎么简化、刚度怎么给、代码怎么写、故障特征怎么加、模型怎么验证。适合做传动系统设计、状态监测算法和故障诊断的工程师与研究生。2. 齿轮动力学模型的核心集中质量法与时变啮合刚度2.1 为什么把齿轮副简化成两个圆盘加一根时变弹簧完整齿轮动力学涉及轮齿弯曲变形、轮体弹性、轴承油膜、箱体结构建立全柔性体模型不仅代价高而且大量参数在工程现场根本测量不到。常见做法是采用集中质量法把一对啮合齿轮简化成两个刚性圆盘由一根刚度随时间变化的弹簧和阻尼器连接。这样系统的自由度从几十万降到一个或几个却仍然保留了齿轮动力学最核心的可观测特征——啮合频率及其谐波激励下的振动响应。对于直齿圆柱齿轮副取啮合线方向的相对位移作为广义坐标运动方程写成m_e · x(t) c_m · x(t) k_m(t) · f(x) F_m F_e(t)m_e 是等效质量m_e (J1·J2) / (J1·r_b2² J2·r_b1²)J 为转动惯量r_b 为基圆半径c_m 是啮合阻尼通常按 c_m 2ζ·sqrt(k_m·m_e) 估算ζ 取 0.030.1k_m(t) 就是时变啮合刚度这是齿轮动力学区别于普通转子动力学的关键f(x) 是间隙函数线性模型里 f(x)x含齿侧间隙时变成分段函数F_m 是平均传递载荷F_e(t) 是传动误差产生的位移激励这个方程虽然只有一个自由度却完整抓住了齿轮振动最主要的激励机制。集中质量法在工程上被验证了几十年齿轮箱故障诊断里常用的边频带分析、共振解调分析其理论基础就是这个模型。自由度太少解决不了的问题比如行星轮均载、轴弯扭耦合需要在后续扩展成多自由度系统但建模思路完全一样。2.2 时变啮合刚度单双齿交替是振动的根源齿轮啮合时啮合点沿啮合线移动参与啮合的轮齿对数在单齿和双齿之间切换。重合度在 12 之间时单个轮齿进入啮合区时是双齿分担载荷到节点附近变成单齿承担全部载荷退出时又恢复双齿。载荷分担比例的周期性变化反映到刚度上就是 k_m(t) 随时间周期性波动。这个波动就是齿轮振动最直接的激励源频率等于轴频乘以齿数也就是啮合频率f_z n · z / 60其中 n 是转速r/minz 是齿数。比如 1500 r/min 下驱动轮 20 齿啮合频率正好是 500 Hz。振动能量主要集中在 f_z 及其 2 倍、3 倍谐波上。刚度的数值怎么给三种常见做法经验公式。ISO 6336 里给出的单齿啮合刚度参考值单位齿宽刚度约 1.1×10⁴1.5×10⁴ N/mm·mm在此基础上乘以齿宽和重合度修正系数。石川公式。把轮齿简化为梯形截面悬臂梁考虑弯曲、剪切和基础弹性变形计算精度高于 ISO 查表法。有限元提取。在 ANSYS 里让一对齿按啮合过程滚动提取接触力与变形量比值。这种方法最准但齿面网格和接触设置耗时适合验证阶段。工程上做快速动力学分析时我一般用傅里叶级数近似k_m(t) k_mean · [1 k_var · cos(2π·f_z·t φ)]k_mean 是平均刚度k_var 是波动率通常 0.150.3重合度越大波动越小φ 是初始相位。更精细的模型可以取前 3 阶谐波用实测刚度曲线拟合系数。建模时不要纠结于某一瞬时的刚度绝对值重要的是波动频率和波动幅值是否合理因为这两个参数直接决定了振动响应的主频率和幅值水平。2.3 阻尼、传动误差与系统固有特性阻尼的取值在齿轮动力学里远比刚度和质量更难确定。啮合阻尼来自油膜挤压和材料内摩擦轴承阻尼和箱体连接阻尼通常需要合并到等效阻尼里。工程上常见做法是先按 0.05 的阻尼比起步再用实测振动数据反推修正。这个参数的敏感性很高阻尼比差一倍共振峰值能差出近一倍所以做完模型后要先做阻尼敏感性分析。传动误差TE定义为实际啮合位置与理论位置的偏差来源于齿形修形、制造公差、装配偏心和轮齿变形。它在方程里表现为右端激励项 F_e(t) k_m(t)·e(t)其中 e(t) 是随转角变化的误差函数。最简单的假设是取为啮合频率下的正弦波e(t) e_amp · sin(2π·f_z·t)幅值从齿形公差标准里查。加上误差激励后系统的响应会同时包含刚度波动和误差波动两类成分的贡献更接近真实工况。做完上述设定后系统的固有频率为 f_n sqrt(k_m / m_e) / (2π)。因为刚度随时间波动严格说这是一个参变系统不是定常系统固有频率也会在小范围内周期性变化。当啮合频率 f_z 接近 f_n 时系统发生主共振振动幅值显著放大f_z 接近 f_n/2 时可能出现次谐波共振。设计阶段的核心任务就是避开这个转速区间这也正是模型要用在的地方。3. 用 MATLAB 建立和求解齿轮动力学模型3.1 把二阶方程改写为状态空间形式MATLAB 的 ode45 只能处理一阶微分方程组所以先把二阶方程改写成状态空间形式。取状态向量 y [x; x]其中 x 是相对位移x 是相对速度得到y1 y2y2 (F_m F_e(t) - c_m·y2 - k_m(t)·f(y1)) / m_e这就是标准的常微分方程初值问题。把所有参数集中到一个参数结构体里用匿名函数或独立函数文件描述右端项。相比把所有变量塞进全局变量用结构体传参的方式在参数扫描时更灵活改一处就能跑一组新工况。初始条件的设置有一个容易被忽略的问题直接从 x(0)0、v(0)0 开始积分系统会经历一段明显的瞬态过程需要几个周期才能进入稳态。如果只关心稳态响应可以把仿真时长拉长到 200 个啮合周期以上然后丢弃前 20% 的数据如果关心的是起动过程或载荷突变响应则要保留完整的瞬态段并单独分析。3.2 ode45 求解的最小可运行代码下面是一段完整可运行的 MATLAB 代码实现单自由度直齿轮副动力学求解。代码里把刚度波动、传动误差和阻尼全部放在一起输出时域响应并做 FFT 频谱分析。% 齿轮动力学建模 -- 单自由度扭转振动模型 % 状态 y [x; x_dot]; x为啮合线相对位移(m), x_dot为相对速度(m/s) clear; clc; close all; % ----- 几何与工况参数 ----- z1 20; % 小齿轮齿数 z2 60; % 大齿轮齿数 n1 1500; % 小齿轮转速(r/min) m 3e-3; % 模数(m), 3mm alpha 20*pi/180; % 压力角(rad) % ----- 等效质量计算 ----- J1 0.5 * pi/4 * (m*z1)^2 * 1e3 * (0.1)^2; % 小齿轮转动惯量估算, 齿宽0.1m, 密度近似 J2 0.5 * pi/4 * (m*z2)^2 * 1e3 * (0.1)^2; rb1 m*z1/2 * cos(alpha); rb2 m*z2/2 * cos(alpha); me (J1*J2) / (J1*rb2^2 J2*rb1^2); % 等效质量 % ----- 刚度与阻尼 ----- k_mean 1.2e8; % 平均啮合刚度(N/m), 按ISO经验估算 k_var 0.20; % 刚度波动率, 重合度1.6~1.8经验值 zeta 0.05; % 阻尼比 cm 2*zeta*sqrt(k_mean*me); % ----- 激励参数 ----- fz n1/60 * z1; % 啮合频率(Hz) load_amp 500; % 单位齿宽载荷(N/mm)乘齿宽后的平均载荷(N) Fm load_amp * 0.1; % 平均传递载荷 e_amp 20e-6; % 传动误差幅值(20um) % ----- 时间轴与初始条件 ----- fs 20000; % 采样频率(Hz), 足够覆盖5倍啮合频率 t_end 0.2; % 仿真时长(s), 约25个轴周期 tspan (0:1/fs:t_end); y0 [0; 0]; % 从静止开始 % ----- 定义ODE右端项 ----- odefun (t, y) gear_ode(t, y, me, cm, k_mean, k_var, fz, e_amp, Fm); % ----- 求解 ----- opts odeset(RelTol, 1e-6, AbsTol, 1e-8, MaxStep, 1/fz/20); [t, y] ode45(odefun, tspan, y0, opts); % ----- 结果可视化 ----- figure; plot(t, y(:,1)*1e6); % 位移单位换算为um xlabel(时间 (s)); ylabel(相对位移 (um)); title(齿轮副啮合线相对位移响应); figure; N length(t); Y fft(y(:,1)); f_ax (0:N-1)/N * fs; amp abs(Y) / N * 2; semilogy(f_ax(1:N/2), amp(1:N/2)); xlim([0 5000]); xlabel(频率 (Hz)); ylabel(幅值 (m)); title(位移响应频谱);对应的子函数gear_ode.mfunction dydt gear_ode(t, y, me, cm, k_mean, k_var, fz, e_amp, Fm) % 时变啮合刚度, 一阶傅里叶近似 kt k_mean * (1 k_var * cos(2*pi*fz*t)); % 传动误差位移激励 et e_amp * sin(2*pi*fz*t); % 激励力 平均载荷 误差引起的动态附加力 F_total Fm kt * et; % 状态方程 dydt zeros(2,1); dydt(1) y(2); dydt(2) (F_total - cm*y(2) - kt*y(1)) / me; end代码逻辑说明主脚本把几何参数折算成等效质量me啮合频率由转速乘齿数得到。ode45 在每个积分步调用gear_ode该函数内部实时计算该时刻的啮合刚度和传动误差组装成当前激励力返回状态导数。MaxStep被限制为啮合周期的 1/20保证刚度波动在积分步内被充分采样如果去掉这个限制ode45 的变步长可能跨过刚度的快速变化段导致高频成分丢失。参数说明k_var是最敏感的调节旋钮。把它设成 0.15 对应高重合度齿轮的平稳啮合设成 0.3 对应低重合度或齿廓偏差较大的齿轮响应幅值会明显上升。e_amp的单位是米制造精度从 7 级到 5 级齿形误差大约从 20μm 到 8μm修改时注意量纲统一。3.3 时域到频域正确读取啮合频率及其边带从仿真结果做 FFT 时三个问题最常见。第一频率分辨率。分辨率等于 fs/N仿真时长越短谱线越粗啮合频率附近的两条边带可能被糊成一条。建议分辨率至少达到 1 Hz也就是说 0.2s 时长对应 5 Hz 分辨率如果要做边频带分析需要把 t_end 加到 1s 以上。第二泄漏。直接 FFT 矩形窗会引入严重的频谱泄漏振动分析至少用 Hanning 窗做边带分析时用 Flat Top 窗更合适。对上述代码中的amp计算处加window hanning(N)然后对窗函数做幅值修正。第三坐标单位。时域位移是米频谱幅值直接输出也是米现场习惯看 mm/s 或 m/s²需要先对位移数值微分一次得到速度、两次得到加速度再做 FFT。MATLAB 里用diff(y(:,1))*fs做粗微分会引入噪声推荐在时域先用sgolay滤波器平滑再微分。啮合频率处的峰值几乎必然存在真正有诊断价值的是它的边带结构。如果在大齿轮转频 f_r2 处出现调制边带也就是 f_z ± k·f_r2说明该齿轮存在安装偏心或齿距累积误差边带间隔越小故障源越可能是低速轴上的局部损伤。这种情况下把频谱横轴放大到 f_z 附近 ±200 Hz 范围观察。4. 从模型到工程结论齿侧间隙、局部故障与参数敏感性4.1 齿侧间隙的分段函数建模与冲击响应实际齿轮副为保证润滑和热膨胀预留了侧隙侧隙让啮合力不能全程保持线性传递。齿轮在换向时齿面先脱离接触再撞向另一侧齿面产生周期性冲击。在线性模型里把这个现象描述的替换方案是把刚度-位移关系改成三段式x b齿面接触恢复力 k(t)·(x - b)-b ≤ x ≤ b脱齿恢复力为 0x -b反向接触恢复力 k(t)·(x b)function F_spring backlash_force(x, b, kt) % x: 相对位移(m), b: 齿侧间隙半宽(m), kt: 当前时刻刚度 if x b F_spring kt * (x - b); % 正向接触 elseif x -b F_spring kt * (x b); % 反向接触 else F_spring 0; % 空行程 end end这段函数把原方程里的kt * y(1)替换掉即可。注意b是半间隙齿轮手册里标注的侧隙通常指圆周方向的游隙沿啮合线方向的当量间隙需要除以压力角的余弦。加入间隙后在时域响应里会观察到明显的“敲击”脉冲频谱上的表现是宽频背景噪声升高啮合频率处的峰值反而不如线性模型突出。这时候如果还按线性模型的思路去查单一特征频率容易漏判。间隙系统的另一个特征是跳跃现象升速和降速过程中在临界转速附近振动幅值不沿同一路径变化出现迟滞。做转速扫描时要注意区分升速和降速两组曲线。4.2 局部损伤的刚度损失建模裂纹与剥落的边频带齿面剥落和齿根裂纹都会造成局部刚度损失而且这种损失是周期性的——受损轮齿每转入啮合区一次刚度就跌一次。建模时分三步确定受损齿所在轴的转频 f_r也就是 n/60。构造一个刚度缺口函数 g(t)每个轴周期内下降一次缺口宽度对应齿面上损伤区的角度范围通常取啮合周期的 10%30%。将总刚度调整为 k(t) k_healthy(t) - Δk·g(t)其中 Δk 是损伤导致的刚度下降幅值。% 构造局部刚度缺口, 每个轴周期出现一次 fr n1/60; % 小齿轮转频 gap_width 0.2/fz; % 缺口占啮合周期的20% t_local mod(t, 1/fr); % 当前轴周期内的时间 if t_local gap_width g 1 - 0.3 * sin(pi * t_local / gap_width); % 30%刚度损失, 平滑过渡 else g 1; end kt k_mean * (1 k_var * cos(2*pi*fz*t)) * g;这个模型的输出频谱上会出现以故障齿轮转频为间隔的边频带边带数量一般 35 对幅值包络的形状和缺口深度相关。局部故障初期的特征是边带幅值低、靠近中心频率随着损伤扩展边带幅值上升且向远离中心频率的方向散开。做在线监测时可以计算边带与中心频率的幅值比作为趋势指标。这里的建模假设是损伤只改变刚度、不改变载荷实际磨损初期成立严重剥落时需要考虑载荷冲击对相邻齿的影响。4.3 参数敏感性转速扫描与共振区定位模型建立后第一件该做的事不是追求高精度而是跑参数扫描摸清系统在什么转速下会共振。下面这段代码对转速做扫描记录每个转速下的稳态 RMS 振动速度直接输出共振曲线% 转速扫描: 找出共振转速区 n_range 600:50:2400; % 扫描转速范围(r/min) rms_result zeros(size(n_range)); for i 1:length(n_range) n1 n_range(i); fz_i n1/60 * z1; % 重新构造ODE并求解, 保持激励力和阻尼不变 odefun_i (t,y) gear_ode(t, y, me, cm, k_mean, k_var, fz_i, e_amp, Fm); [t_i, y_i] ode45(odefun_i, tspan, y0, opts); % 丢弃瞬态段, 取最后100个啮合周期计算速度RMS n_discard round(length(t_i)*0.2); v_ss diff(y_i(n_discard:end,1)) * fs; rms_result(i) sqrt(mean(v_ss.^2)); end plot(n_range, rms_result*1000, o-); xlabel(转速 (r/min)); ylabel(振动速度 RMS (mm/s));扫描结果通常呈现一个或两个明显的峰峰值转速对应啮合频率接近系统固有频率。如果共振区落在常用工作转速内修改策略按优先级排序调整阻尼比增大ζ实现结构阻尼或改用低阻尼材料则相反、改变重合度通过变位系数调整、改变齿数或模数从而改变 f_z 避开 f_n。这套流程在齿轮设计阶段跑一遍比样机测试后再补救成本低得多。下面是一组典型参数的敏感性规律供调整方向时参考参数调大后的影响调小后的影响工程手段重合度 ε刚度波动率下降响应幅值降低波动率上升啮合冲击增强变位系数、齿顶高系数调整阻尼比 ζ共振峰削平但传递到箱体的高频力衰减有限峰变尖稳态幅值升高结构阻尼、混合材料、油膜特性模数 m刚度增大固有频率上升刚度下降齿根应力增大重新选型时权衡传动误差 e_amp误差激励线性放大幅值上升幅值下降但修形过度导致接触区变窄齿廓修形量、磨齿精度齿侧间隙 b冲击增强脱齿区间变宽噪声变大可能造成干涉和润滑膜破裂装配工艺与热膨胀补偿表格里的规律都对应着明确的物理机制重合度影响激励的“平滑度”阻尼影响共振峰的形态误差是线性输入的激励源。做工程设计时参数调整的优先级应当先看激励源重合度和误差再看系统阻尼最后才动质量和刚度。5. 验证模型可信度的四个手段模型搭建完成后先做验证再做参数扫描顺序不能反过来。验证思路不是拿仿真曲线和现场测点对比那是标定验证要先证明代码本身的数学实现没有错误。第一性验证是去掉所有外激励和阻尼只保留刚度项给系统一个初始位移后让它自由振动。用k_mean和me手算出理论固有频率 f_n sqrt(k_mean/me)/(2π)再对自由衰减响应的位移做 FFT看谱峰是否落在理论值上误差超过 0.5% 说明状态方程组装有问题。第二性验证是能量检查稳态段激励力做的功应该等于阻尼耗散的能量用trapz对功率做时间积分两者相对偏差小于 5% 即合理偏差过大通常意味着时间步长对刚度快速波动阶段的采样不足此时调小MaxStep而不是加大RelTol。第三性验证是收敛性检查把RelTol从 1e-3 往下缩到 1e-6每次计算稳态 RMS如果相邻两档结果变化不到 1%认为当前容差不主导误差。第四性验证是残差回代用deval设立插值后的状态和时间代回odefun计算残差残差量级应远小于加速度特征量这一步能抓出状态方程里量纲不匹配的隐性错误。一个实操技巧把上述四项验证写成一个validate_model.m脚本每次修改参数后先跑验证再跑扫描。模型可信度从“ode45 解出来了”提升到“解对了”这个距离往往就是工程判断是否可靠的差距。收敛性检查尤其值得保留在脚本里因为 MATLAB 不同版本之间求解器执行细节有小差异换环境后模型结果可能变化验证脚本能立刻暴露问题。本文还有配套的精品资源点击获取
基于MATLAB的齿轮动力学建模与振动响应分析
简介面向机械工程与故障诊断领域的MATLAB源码包聚焦齿轮动力学建模与振动信号分析适合需要借助短时傅里叶变换STFT研究齿轮啮合特性、识别齿面磨损或裂纹等故障模式的研究人员和工程师。压缩包共三个文件包括两个.m脚本与一个.asv自动备份文件整体仅5KB体量精巧便于直接阅读和修改。已有3137人学习下载。内容围绕齿轮动力学模型展开通过MATLAB信号处理工具箱中的spectrogram函数实现加窗、切割与可视化可帮助用户掌握从齿轮振动数据中提取时频特征、定位瞬态故障频率的分析思路同时结合窗函数选择、频率分辨率权衡等细节为后续扩展小波变换或有限元分析提供了基础。这套源码对于入门齿轮动力学数值实验和快速搭建振动信号处理流程具有实用参考价值。1. 齿轮动力学建模从一对齿的啮合到整条传动链的振动响应齿轮箱振动超标、减速机异响、齿面早期点蚀这些现场问题背后几乎都指向同一个源头——啮合副的动态激励。手算齿轮强度只能给出静载安全系数回答不了“为什么这台箱体在某个转速下噪声特别大”这类问题。要回答就得建立齿轮动力学模型把齿轮副的质量、刚度、阻尼和激励放到一个微分方程组里用数值方法求解时域响应再从响应里找共振点、边频带和冲击特征。MATLAB在这条路上几乎是默认工具ode45做数值积分、FFT做频谱分析、优化工具箱做参数识别一个环境全链打通。这篇文章按常见做法讲一套能落地的扭转振动模型理论怎么简化、刚度怎么给、代码怎么写、故障特征怎么加、模型怎么验证。适合做传动系统设计、状态监测算法和故障诊断的工程师与研究生。2. 齿轮动力学模型的核心集中质量法与时变啮合刚度2.1 为什么把齿轮副简化成两个圆盘加一根时变弹簧完整齿轮动力学涉及轮齿弯曲变形、轮体弹性、轴承油膜、箱体结构建立全柔性体模型不仅代价高而且大量参数在工程现场根本测量不到。常见做法是采用集中质量法把一对啮合齿轮简化成两个刚性圆盘由一根刚度随时间变化的弹簧和阻尼器连接。这样系统的自由度从几十万降到一个或几个却仍然保留了齿轮动力学最核心的可观测特征——啮合频率及其谐波激励下的振动响应。对于直齿圆柱齿轮副取啮合线方向的相对位移作为广义坐标运动方程写成m_e · x(t) c_m · x(t) k_m(t) · f(x) F_m F_e(t)m_e 是等效质量m_e (J1·J2) / (J1·r_b2² J2·r_b1²)J 为转动惯量r_b 为基圆半径c_m 是啮合阻尼通常按 c_m 2ζ·sqrt(k_m·m_e) 估算ζ 取 0.030.1k_m(t) 就是时变啮合刚度这是齿轮动力学区别于普通转子动力学的关键f(x) 是间隙函数线性模型里 f(x)x含齿侧间隙时变成分段函数F_m 是平均传递载荷F_e(t) 是传动误差产生的位移激励这个方程虽然只有一个自由度却完整抓住了齿轮振动最主要的激励机制。集中质量法在工程上被验证了几十年齿轮箱故障诊断里常用的边频带分析、共振解调分析其理论基础就是这个模型。自由度太少解决不了的问题比如行星轮均载、轴弯扭耦合需要在后续扩展成多自由度系统但建模思路完全一样。2.2 时变啮合刚度单双齿交替是振动的根源齿轮啮合时啮合点沿啮合线移动参与啮合的轮齿对数在单齿和双齿之间切换。重合度在 12 之间时单个轮齿进入啮合区时是双齿分担载荷到节点附近变成单齿承担全部载荷退出时又恢复双齿。载荷分担比例的周期性变化反映到刚度上就是 k_m(t) 随时间周期性波动。这个波动就是齿轮振动最直接的激励源频率等于轴频乘以齿数也就是啮合频率f_z n · z / 60其中 n 是转速r/minz 是齿数。比如 1500 r/min 下驱动轮 20 齿啮合频率正好是 500 Hz。振动能量主要集中在 f_z 及其 2 倍、3 倍谐波上。刚度的数值怎么给三种常见做法经验公式。ISO 6336 里给出的单齿啮合刚度参考值单位齿宽刚度约 1.1×10⁴1.5×10⁴ N/mm·mm在此基础上乘以齿宽和重合度修正系数。石川公式。把轮齿简化为梯形截面悬臂梁考虑弯曲、剪切和基础弹性变形计算精度高于 ISO 查表法。有限元提取。在 ANSYS 里让一对齿按啮合过程滚动提取接触力与变形量比值。这种方法最准但齿面网格和接触设置耗时适合验证阶段。工程上做快速动力学分析时我一般用傅里叶级数近似k_m(t) k_mean · [1 k_var · cos(2π·f_z·t φ)]k_mean 是平均刚度k_var 是波动率通常 0.150.3重合度越大波动越小φ 是初始相位。更精细的模型可以取前 3 阶谐波用实测刚度曲线拟合系数。建模时不要纠结于某一瞬时的刚度绝对值重要的是波动频率和波动幅值是否合理因为这两个参数直接决定了振动响应的主频率和幅值水平。2.3 阻尼、传动误差与系统固有特性阻尼的取值在齿轮动力学里远比刚度和质量更难确定。啮合阻尼来自油膜挤压和材料内摩擦轴承阻尼和箱体连接阻尼通常需要合并到等效阻尼里。工程上常见做法是先按 0.05 的阻尼比起步再用实测振动数据反推修正。这个参数的敏感性很高阻尼比差一倍共振峰值能差出近一倍所以做完模型后要先做阻尼敏感性分析。传动误差TE定义为实际啮合位置与理论位置的偏差来源于齿形修形、制造公差、装配偏心和轮齿变形。它在方程里表现为右端激励项 F_e(t) k_m(t)·e(t)其中 e(t) 是随转角变化的误差函数。最简单的假设是取为啮合频率下的正弦波e(t) e_amp · sin(2π·f_z·t)幅值从齿形公差标准里查。加上误差激励后系统的响应会同时包含刚度波动和误差波动两类成分的贡献更接近真实工况。做完上述设定后系统的固有频率为 f_n sqrt(k_m / m_e) / (2π)。因为刚度随时间波动严格说这是一个参变系统不是定常系统固有频率也会在小范围内周期性变化。当啮合频率 f_z 接近 f_n 时系统发生主共振振动幅值显著放大f_z 接近 f_n/2 时可能出现次谐波共振。设计阶段的核心任务就是避开这个转速区间这也正是模型要用在的地方。3. 用 MATLAB 建立和求解齿轮动力学模型3.1 把二阶方程改写为状态空间形式MATLAB 的 ode45 只能处理一阶微分方程组所以先把二阶方程改写成状态空间形式。取状态向量 y [x; x]其中 x 是相对位移x 是相对速度得到y1 y2y2 (F_m F_e(t) - c_m·y2 - k_m(t)·f(y1)) / m_e这就是标准的常微分方程初值问题。把所有参数集中到一个参数结构体里用匿名函数或独立函数文件描述右端项。相比把所有变量塞进全局变量用结构体传参的方式在参数扫描时更灵活改一处就能跑一组新工况。初始条件的设置有一个容易被忽略的问题直接从 x(0)0、v(0)0 开始积分系统会经历一段明显的瞬态过程需要几个周期才能进入稳态。如果只关心稳态响应可以把仿真时长拉长到 200 个啮合周期以上然后丢弃前 20% 的数据如果关心的是起动过程或载荷突变响应则要保留完整的瞬态段并单独分析。3.2 ode45 求解的最小可运行代码下面是一段完整可运行的 MATLAB 代码实现单自由度直齿轮副动力学求解。代码里把刚度波动、传动误差和阻尼全部放在一起输出时域响应并做 FFT 频谱分析。% 齿轮动力学建模 -- 单自由度扭转振动模型 % 状态 y [x; x_dot]; x为啮合线相对位移(m), x_dot为相对速度(m/s) clear; clc; close all; % ----- 几何与工况参数 ----- z1 20; % 小齿轮齿数 z2 60; % 大齿轮齿数 n1 1500; % 小齿轮转速(r/min) m 3e-3; % 模数(m), 3mm alpha 20*pi/180; % 压力角(rad) % ----- 等效质量计算 ----- J1 0.5 * pi/4 * (m*z1)^2 * 1e3 * (0.1)^2; % 小齿轮转动惯量估算, 齿宽0.1m, 密度近似 J2 0.5 * pi/4 * (m*z2)^2 * 1e3 * (0.1)^2; rb1 m*z1/2 * cos(alpha); rb2 m*z2/2 * cos(alpha); me (J1*J2) / (J1*rb2^2 J2*rb1^2); % 等效质量 % ----- 刚度与阻尼 ----- k_mean 1.2e8; % 平均啮合刚度(N/m), 按ISO经验估算 k_var 0.20; % 刚度波动率, 重合度1.6~1.8经验值 zeta 0.05; % 阻尼比 cm 2*zeta*sqrt(k_mean*me); % ----- 激励参数 ----- fz n1/60 * z1; % 啮合频率(Hz) load_amp 500; % 单位齿宽载荷(N/mm)乘齿宽后的平均载荷(N) Fm load_amp * 0.1; % 平均传递载荷 e_amp 20e-6; % 传动误差幅值(20um) % ----- 时间轴与初始条件 ----- fs 20000; % 采样频率(Hz), 足够覆盖5倍啮合频率 t_end 0.2; % 仿真时长(s), 约25个轴周期 tspan (0:1/fs:t_end); y0 [0; 0]; % 从静止开始 % ----- 定义ODE右端项 ----- odefun (t, y) gear_ode(t, y, me, cm, k_mean, k_var, fz, e_amp, Fm); % ----- 求解 ----- opts odeset(RelTol, 1e-6, AbsTol, 1e-8, MaxStep, 1/fz/20); [t, y] ode45(odefun, tspan, y0, opts); % ----- 结果可视化 ----- figure; plot(t, y(:,1)*1e6); % 位移单位换算为um xlabel(时间 (s)); ylabel(相对位移 (um)); title(齿轮副啮合线相对位移响应); figure; N length(t); Y fft(y(:,1)); f_ax (0:N-1)/N * fs; amp abs(Y) / N * 2; semilogy(f_ax(1:N/2), amp(1:N/2)); xlim([0 5000]); xlabel(频率 (Hz)); ylabel(幅值 (m)); title(位移响应频谱);对应的子函数gear_ode.mfunction dydt gear_ode(t, y, me, cm, k_mean, k_var, fz, e_amp, Fm) % 时变啮合刚度, 一阶傅里叶近似 kt k_mean * (1 k_var * cos(2*pi*fz*t)); % 传动误差位移激励 et e_amp * sin(2*pi*fz*t); % 激励力 平均载荷 误差引起的动态附加力 F_total Fm kt * et; % 状态方程 dydt zeros(2,1); dydt(1) y(2); dydt(2) (F_total - cm*y(2) - kt*y(1)) / me; end代码逻辑说明主脚本把几何参数折算成等效质量me啮合频率由转速乘齿数得到。ode45 在每个积分步调用gear_ode该函数内部实时计算该时刻的啮合刚度和传动误差组装成当前激励力返回状态导数。MaxStep被限制为啮合周期的 1/20保证刚度波动在积分步内被充分采样如果去掉这个限制ode45 的变步长可能跨过刚度的快速变化段导致高频成分丢失。参数说明k_var是最敏感的调节旋钮。把它设成 0.15 对应高重合度齿轮的平稳啮合设成 0.3 对应低重合度或齿廓偏差较大的齿轮响应幅值会明显上升。e_amp的单位是米制造精度从 7 级到 5 级齿形误差大约从 20μm 到 8μm修改时注意量纲统一。3.3 时域到频域正确读取啮合频率及其边带从仿真结果做 FFT 时三个问题最常见。第一频率分辨率。分辨率等于 fs/N仿真时长越短谱线越粗啮合频率附近的两条边带可能被糊成一条。建议分辨率至少达到 1 Hz也就是说 0.2s 时长对应 5 Hz 分辨率如果要做边频带分析需要把 t_end 加到 1s 以上。第二泄漏。直接 FFT 矩形窗会引入严重的频谱泄漏振动分析至少用 Hanning 窗做边带分析时用 Flat Top 窗更合适。对上述代码中的amp计算处加window hanning(N)然后对窗函数做幅值修正。第三坐标单位。时域位移是米频谱幅值直接输出也是米现场习惯看 mm/s 或 m/s²需要先对位移数值微分一次得到速度、两次得到加速度再做 FFT。MATLAB 里用diff(y(:,1))*fs做粗微分会引入噪声推荐在时域先用sgolay滤波器平滑再微分。啮合频率处的峰值几乎必然存在真正有诊断价值的是它的边带结构。如果在大齿轮转频 f_r2 处出现调制边带也就是 f_z ± k·f_r2说明该齿轮存在安装偏心或齿距累积误差边带间隔越小故障源越可能是低速轴上的局部损伤。这种情况下把频谱横轴放大到 f_z 附近 ±200 Hz 范围观察。4. 从模型到工程结论齿侧间隙、局部故障与参数敏感性4.1 齿侧间隙的分段函数建模与冲击响应实际齿轮副为保证润滑和热膨胀预留了侧隙侧隙让啮合力不能全程保持线性传递。齿轮在换向时齿面先脱离接触再撞向另一侧齿面产生周期性冲击。在线性模型里把这个现象描述的替换方案是把刚度-位移关系改成三段式x b齿面接触恢复力 k(t)·(x - b)-b ≤ x ≤ b脱齿恢复力为 0x -b反向接触恢复力 k(t)·(x b)function F_spring backlash_force(x, b, kt) % x: 相对位移(m), b: 齿侧间隙半宽(m), kt: 当前时刻刚度 if x b F_spring kt * (x - b); % 正向接触 elseif x -b F_spring kt * (x b); % 反向接触 else F_spring 0; % 空行程 end end这段函数把原方程里的kt * y(1)替换掉即可。注意b是半间隙齿轮手册里标注的侧隙通常指圆周方向的游隙沿啮合线方向的当量间隙需要除以压力角的余弦。加入间隙后在时域响应里会观察到明显的“敲击”脉冲频谱上的表现是宽频背景噪声升高啮合频率处的峰值反而不如线性模型突出。这时候如果还按线性模型的思路去查单一特征频率容易漏判。间隙系统的另一个特征是跳跃现象升速和降速过程中在临界转速附近振动幅值不沿同一路径变化出现迟滞。做转速扫描时要注意区分升速和降速两组曲线。4.2 局部损伤的刚度损失建模裂纹与剥落的边频带齿面剥落和齿根裂纹都会造成局部刚度损失而且这种损失是周期性的——受损轮齿每转入啮合区一次刚度就跌一次。建模时分三步确定受损齿所在轴的转频 f_r也就是 n/60。构造一个刚度缺口函数 g(t)每个轴周期内下降一次缺口宽度对应齿面上损伤区的角度范围通常取啮合周期的 10%30%。将总刚度调整为 k(t) k_healthy(t) - Δk·g(t)其中 Δk 是损伤导致的刚度下降幅值。% 构造局部刚度缺口, 每个轴周期出现一次 fr n1/60; % 小齿轮转频 gap_width 0.2/fz; % 缺口占啮合周期的20% t_local mod(t, 1/fr); % 当前轴周期内的时间 if t_local gap_width g 1 - 0.3 * sin(pi * t_local / gap_width); % 30%刚度损失, 平滑过渡 else g 1; end kt k_mean * (1 k_var * cos(2*pi*fz*t)) * g;这个模型的输出频谱上会出现以故障齿轮转频为间隔的边频带边带数量一般 35 对幅值包络的形状和缺口深度相关。局部故障初期的特征是边带幅值低、靠近中心频率随着损伤扩展边带幅值上升且向远离中心频率的方向散开。做在线监测时可以计算边带与中心频率的幅值比作为趋势指标。这里的建模假设是损伤只改变刚度、不改变载荷实际磨损初期成立严重剥落时需要考虑载荷冲击对相邻齿的影响。4.3 参数敏感性转速扫描与共振区定位模型建立后第一件该做的事不是追求高精度而是跑参数扫描摸清系统在什么转速下会共振。下面这段代码对转速做扫描记录每个转速下的稳态 RMS 振动速度直接输出共振曲线% 转速扫描: 找出共振转速区 n_range 600:50:2400; % 扫描转速范围(r/min) rms_result zeros(size(n_range)); for i 1:length(n_range) n1 n_range(i); fz_i n1/60 * z1; % 重新构造ODE并求解, 保持激励力和阻尼不变 odefun_i (t,y) gear_ode(t, y, me, cm, k_mean, k_var, fz_i, e_amp, Fm); [t_i, y_i] ode45(odefun_i, tspan, y0, opts); % 丢弃瞬态段, 取最后100个啮合周期计算速度RMS n_discard round(length(t_i)*0.2); v_ss diff(y_i(n_discard:end,1)) * fs; rms_result(i) sqrt(mean(v_ss.^2)); end plot(n_range, rms_result*1000, o-); xlabel(转速 (r/min)); ylabel(振动速度 RMS (mm/s));扫描结果通常呈现一个或两个明显的峰峰值转速对应啮合频率接近系统固有频率。如果共振区落在常用工作转速内修改策略按优先级排序调整阻尼比增大ζ实现结构阻尼或改用低阻尼材料则相反、改变重合度通过变位系数调整、改变齿数或模数从而改变 f_z 避开 f_n。这套流程在齿轮设计阶段跑一遍比样机测试后再补救成本低得多。下面是一组典型参数的敏感性规律供调整方向时参考参数调大后的影响调小后的影响工程手段重合度 ε刚度波动率下降响应幅值降低波动率上升啮合冲击增强变位系数、齿顶高系数调整阻尼比 ζ共振峰削平但传递到箱体的高频力衰减有限峰变尖稳态幅值升高结构阻尼、混合材料、油膜特性模数 m刚度增大固有频率上升刚度下降齿根应力增大重新选型时权衡传动误差 e_amp误差激励线性放大幅值上升幅值下降但修形过度导致接触区变窄齿廓修形量、磨齿精度齿侧间隙 b冲击增强脱齿区间变宽噪声变大可能造成干涉和润滑膜破裂装配工艺与热膨胀补偿表格里的规律都对应着明确的物理机制重合度影响激励的“平滑度”阻尼影响共振峰的形态误差是线性输入的激励源。做工程设计时参数调整的优先级应当先看激励源重合度和误差再看系统阻尼最后才动质量和刚度。5. 验证模型可信度的四个手段模型搭建完成后先做验证再做参数扫描顺序不能反过来。验证思路不是拿仿真曲线和现场测点对比那是标定验证要先证明代码本身的数学实现没有错误。第一性验证是去掉所有外激励和阻尼只保留刚度项给系统一个初始位移后让它自由振动。用k_mean和me手算出理论固有频率 f_n sqrt(k_mean/me)/(2π)再对自由衰减响应的位移做 FFT看谱峰是否落在理论值上误差超过 0.5% 说明状态方程组装有问题。第二性验证是能量检查稳态段激励力做的功应该等于阻尼耗散的能量用trapz对功率做时间积分两者相对偏差小于 5% 即合理偏差过大通常意味着时间步长对刚度快速波动阶段的采样不足此时调小MaxStep而不是加大RelTol。第三性验证是收敛性检查把RelTol从 1e-3 往下缩到 1e-6每次计算稳态 RMS如果相邻两档结果变化不到 1%认为当前容差不主导误差。第四性验证是残差回代用deval设立插值后的状态和时间代回odefun计算残差残差量级应远小于加速度特征量这一步能抓出状态方程里量纲不匹配的隐性错误。一个实操技巧把上述四项验证写成一个validate_model.m脚本每次修改参数后先跑验证再跑扫描。模型可信度从“ode45 解出来了”提升到“解对了”这个距离往往就是工程判断是否可靠的差距。收敛性检查尤其值得保留在脚本里因为 MATLAB 不同版本之间求解器执行细节有小差异换环境后模型结果可能变化验证脚本能立刻暴露问题。本文还有配套的精品资源点击获取