简介面向机械工程与车辆工程领域的齿轮动力学学习者与研究者资源围绕22自由度齿轮传动系统的建模与数值仿真展开重点解决由齿侧间隙、时变啮合刚度等因素引发的非线性振动响应计算问题。压缩包为rar格式大小仅3KB内含1个ToqurVibratory.m脚本文件。脚本基于MATLAB平台利用ode45龙格-库塔积分器对系统状态方程进行求解涵盖了状态变量设定、质量/刚度/阻尼矩阵组装、外载荷定义、初值条件输入以及结果后处理等关键环节并用图形方式展示角位移、角速度等动态曲线。目前已有300人学习浏览。对于正在学习机械振动、齿轮动力学或MATLAB数值分析的人员这份轻量代码提供了可直接加载运行的参考模板便于深入理解非线性微分方程的求解流程并可在其基础上修改参数、扩展自由度快速开展齿轮传动系统的动态特性研究与优化验证。1. 从啮合冲击到 22 自由度ToqurVibratory 模型的切入点一台减速机齿面出现点蚀时现场振动频谱往往要到故障中后期才显露出边带。最麻烦的是你很难区分是齿形误差还是轴弯曲在激励齿轮。要回答这个问题与其对着频谱猜不如把齿轮系统建成多自由度动力学模型用仿真把啮合冲击、时变刚度和轴系扭转放到同一套方程里看。ToqurVibratory.rar 里的 ToqurVibratory.m 正是这样的实现它用 MATLAB 的 ode45 求解 22 自由度齿轮动力学方程输出角位移、角速度等状态量。这个压缩包适合传动设计、故障诊断工程师也适合想把动力学方程真正跑起来的研究生。读懂它你就掌握了一整套从建模、求解到结果的齿轮动力学工作流。2. 22 自由度齿轮模型的坐标系、自由度分配与方程组装2.1 自由度的物理含义与分配齿轮动力学计算的第一步是决定模型中哪些部件能“动”。以多级平行轴齿轮箱为例每一根轴的转动、齿轮副啮合点的相对位移、轴承座的弹性位移甚至电机转子的扭转都可能对动态响应产生影响。22 自由度意味着系统已经做成一个集中质量网络每个广义坐标对应一个弹性势能或动能项而不是把齿轮当成刚体简单折算。在 ToqurVibratory 这类模型里自由度通常按部件分组。表 2-1 是一种常见分配方式具体到某个压缩包会因齿轮级数和轴段划分而异但整体思路一致。表 2-1 22 个自由度的分组示例自由度区间物理量对应部件12角位移、横向位移电机转子 / 输入小齿轮38角位移、水平和垂直位移中间轴大齿轮、小齿轮916角位移、水平/垂直/轴向位移输出轴齿轮副与轴承座1722齿面法向相对位移、传递误差协调量多对啮合齿轮副的啮合点这里需要说明12 只是一个区间示例实际如果输入轴是一根刚性轴可能只有一个扭转自由度。但 22 自由度模型的每一列位移都会进入质量矩阵和刚度矩阵直接决定求解规模。自由度的取舍原则很简单能引起齿面载荷波动的弹性变形必须保留比如轴的扭转变形和齿轮副的啮合变形对结果影响很小的刚体平动可以合并或固定。2.2 二阶微分方程与状态变换齿轮系统的运动方程通常写成M q C q K(t) q T(t)式中 q 是 22 维广义位移向量q 和 q 是对应速度与加速度。M 是质量矩阵C 是阻尼矩阵K(t) 是时变刚度矩阵T(t) 是外载荷。K(t) 随时间变化是因为啮合点位置随齿轮转角移动齿对啮合刚度的变化是齿轮动力学区别于普通转子动力学的重要特征。ode45 要求的是常微分方程组的初值问题所以必须把二阶方程转换成一阶状态方程。常见做法是令 y [q; q]得到 44 维状态向量。下面代码展示如何在 MATLAB 函数里定义这个过程。function dydt gear22(t, y, p) % 22自由度齿轮系统的一阶状态方程 % y(1:22) 为广义位移, y(23:44) 为广义速度 q y(1:22); v y(23:44); % 时变啮合刚度: 平均刚度叠加一阶余弦波动 % p.z 齿数, q(1) 为主动轮角位移 k_m p.k_mean p.k_amp * cos(p.z * q(1) p.phi0); % 刚度矩阵中啮合位置耦合项 K p.K0; K(p.pinion_dof, p.gear_dof) K(p.pinion_dof, p.gear_dof) - k_m; K(p.gear_dof, p.pinion_dof) K(p.gear_dof, p.pinion_dof) - k_m; K(p.pinion_dof, p.pinion_dof) K(p.pinion_dof, p.pinion_dof) k_m; K(p.gear_dof, p.gear_dof) K(p.gear_dof, p.gear_dof) k_m; % 外力: 输入扭矩, 负载扭矩, 啮合误差激励 T zeros(22, 1); T(p.pinion_dof) p.T_in; T(p.gear_dof) -p.T_load; % M q T - C v - K q accel p.M \ (T - p.C * v - K * q); dydt [v; accel]; end这段代码的逻辑是先取出位移和速度再根据当前角度更新啮合刚度。刚度矩阵中同一组啮合自由度会出现 k_m 和 -k_m本质上是把两齿轮通过弹簧连接起来。p.pinion_dof和p.gear_dof是预先存好的自由度编号比如 2 和 3这样写比直接填数字更不容易出错。这里的参数p.k_mean是平均啮合刚度单位 N/mp.k_amp是刚度波动幅值通常取平均值的 10%30%。p.phi0是初始相位用来对齐齿距误差。注意M \是 MATLAB 左除它对稀疏矩阵尤其高效如果 M 是常矩阵可以在主脚本里提前做p.invM inv(M)然后用accel p.invM * (...)但要注意逆矩阵会丢失稀疏性22 自由度规模无所谓更大模型建议保留左除。2.3 参数表与激励设置要真正让模型算得动手里的几何参数必须换算成 M、C、K。表 2-2 给出一组常见示例值便于调试时对照。表 2-2 齿轮动力学计算常用参数示例参数符号示例值说明齿数z25齿轮副主从动轮独立给平均啮合刚度k_mean8e7 N/m按 ISO 6336 查表或有限元计算刚度波动幅值k_amp1.6e7 N/m一般是平均刚度的 20%阻尼比ξ0.03钢-钢齿轮副的经验值输入扭矩T_in50 N·m电机额定输出输入转速n_in1500 rpm决定啮合周期把表里的刚度除以质量矩阵对应元素就能估算系统的固有频率。例如某个齿轮副等效质量 0.05 kg啮合刚度 8e7 N/m单自由度固有频率约为 6360 rad/s约 1012 Hz。这个频率远高于转频但低于采样率所以后处理时要用足够小的最大步长才能捕获。外载荷不一定是常数。实际工况里负载突变会造成扭转冲击ToqurVibratory 这类模型在T(22) -p.T_load的位置可以改为时间分段函数。比如模拟加载瞬间可以用p.T_load * (t 0.1)让负载在 0.1 秒后加入观察瞬态冲击响应。3. ode45 求解 22 自由度系统的步长控制与刚性问题3.1 ode45 的适用边界与方法特点ode45 是 MATLAB 内置的自适应龙格-库塔法它同时使用四阶和五阶公式通过两者差估算局部误差自动调整步长。对齿轮动力学模型啮合刚度随时间连续变化不存在剧烈的尺度分裂所以 ode45 往往比隐式方法更快。但 22 自由度系统的固有频率分散如果刚度阵里同时有 1e4 和 1e8 N/m 的数量级局部误差控制会迫使步长变小仿真时间猛增。这时首先要检查的是单位而不是怀疑求解器。常用做法是先把齿轮副的啮合频率算出来用1/(20*fm)作为最大步长上界。这样不管误差控制怎么变化高频成分至少每个周期采 20 个点。为了保证 FFT 时频率分辨率够用还要保证仿真时间至少 200 个啮合周期。3.2 ode45 的标准调用方式在 MATLAB 主脚本中调用 gear22 的方式如下。这段代码可以直接替换 ToqurVibratory.m 里的求解部分前提是你有对应的 p 参数结构体。% ---- 主求解脚本 ---- p.M 0.05 * eye(22); % 演示用对角质量阵, 实际按转动惯量计算 p.C 0.03 * eye(22); % 比例阻尼, 单位 N·s/m p.K0 5e6 * eye(22); % 基础刚度阵, 单位 N/m p.k_mean 8e7; % 啮合平均刚度 p.k_amp 1.6e7; % 啮合刚度波动 p.z 25; % 主动轮齿数 p.phi0 0; % 初始相位 p.pinion_dof 2; % 主动轮自由度编号 p.gear_dof 3; % 从动轮自由度编号 p.T_in 50; % 输入扭矩 N·m p.T_load 50; % 负载扭矩 N·m y0 zeros(44, 1); % 初始静止 tspan [0 1]; % 仿真 1 秒 opts odeset(RelTol, 1e-7, AbsTol, 1e-9, MaxStep, 5e-5); [t, Y] ode45((t, y) gear22(t, y, p), tspan, y0, opts);这里的关键是AbsTol设置成 1e-9因为齿轮位移量级很小尤其是齿面法向变形可能只有微米级如果只调RelTol绝对误差过大会把微变形直接抹平。MaxStep5e-5在 1500 rpm、齿数 25 的工况下对应啮合频率 625 Hz步长 5e-5 秒时每个啮合周期约 32 步足够看清主要振动。3.3 求解精度与效率的平衡技巧表 3-1 给出 ode45 常用参数推荐这个表也适用于其他刚度和质量量级的齿轮模型。表 3-1 ode45 求解参数推荐参数推荐范围调整依据RelTol1e-6 ~ 1e-8需要捕捉边带时收紧AbsTol1e-8 ~ 1e-10位移量级为微米时取 1e-10MaxStep1/(20fm) ~ 1/(50fm)啮合频率高时取小值InitialStep自动大部分情况不需要手动工程上我的习惯是先用RelTol1e-6跑通再收紧到 1e-8比较两次结果的幅值差。如果幅值差超过 5%说明模型可能进入了刚性问题区域。判断方法很简单把opts.Stats打开看 ode45 是否频繁拒绝步长。如果 rejected 数量达到 accepted 数量的一半以上就该考虑换 ode15s。3.4 刚性问题排查与求解器切换齿轮动力学中的刚性问题通常来自轴承油膜刚度或齿面接触刚度的数值量级差异。比如齿面刚度 1e8 N/m而扭振轴段刚度 1e5 N·m/rad这会让雅可比矩阵特征值相差 3 个数量级。ode45 为了稳定会不断减小步长最终表现为卡在某个时间点不动。遇到这种情况先做三件事。第一把质量矩阵和刚度矩阵各元素的单位统一确认没有 N 和 N·m 混用。第二把阻尼比调大到 0.05 再试增加阻尼可以平滑瞬态但不解决本质问题。第三直接切到ode15s只需要把ode45换成ode15s其余代码不动[t, Y] ode15s((t, y) gear22(t, y, p), tspan, y0, opts);ode15s 是隐式变阶求解器适合刚性系统但它每个步长要解线性方程组22 自由度规模下代价很低。切过去后要重新验证时间序列是否与 ode45 在小步长下收敛到同一曲线。通常两者在驱动转矩恒定且不存在突变时差异很小如果波动很大往往不是求解器的问题而是刚度矩阵里出现了额外自由度冲突比如两个齿轮自由度重复表达了同一个位移。4. 从时间序列到频谱边带ToqurVibratory 结果的后处理流程4.1 变步长数据的等间隔采样ode45 返回的 t 和 Y 是非等间隔的。要计算 FFT必须先把信号重采样到等间隔时间轴。这里我一般选fs 10000Hz能覆盖啮合频率及其高次谐波。用interp1插值时要避免spline产生过冲pchip更适合振动信号。% 取第 8 个自由度的横向振动速度 fs 10000; t_end t(end); tq 0:1/fs:t_end; Vq interp1(t, Y(:, 822), tq, pchip); % 注意索引, 位移22 才是速度 % 去除直流分量, 便于观察振动幅值 Vq Vq - mean(Vq); figure; plot(tq, Vq * 1000); xlabel(时间 (s)); ylabel(振动速度 (mm/s));这里取的是Y(:, 822)因为前 22 列是位移后 22 列是速度。Vq单位是 m/s乘以 1000 得到 mm/s方便与现场测点数据比较。如果看的是角速度单位是 rad/s后处理时不需要换算直接做频谱即可。4.2 频率谱计算与啮合频率定位等间隔采样后用 fft 计算单边幅值谱。下面这段代码会去掉镜像部分并只显示前 3000 Hz 的频率成分一般足够覆盖齿轮故障相关频带。L length(Vq); Yf fft(Vq); P2 abs(Yf / L); P1 P2(1:L/21); P1(2:end-1) 2 * P1(2:end-1); f fs * (0:(L/2)) / L; % 只显示 3000 Hz 以内 idx f 3000; figure; plot(f(idx), P1(idx)); xlabel(频率 (Hz)); ylabel(幅值);注意P1(2:end-1)2*P1(...)是单边谱的标准做法直流分量不需要加倍。仿真数据没有传感器噪声所以频谱会比实测干净很多偶尔会出现非常尖锐的谱线这是正常现象。实际分析时要把频率分辨率dffs/L计算出来如果df大于 5 Hz可能需要延长仿真时间或补零来提高谱线密度。4.3 齿轮振动特征频率与边带识别在齿轮动力学仿真中最常看的三个频率成分是轴转频、啮合频率和边带。表 4-1 列出它们的计算方式。表 4-1 齿轮动力学特征频率计算特征公式示例值(1500rpm, z25)输入轴转频f_r n/6025 Hz啮合频率f_m z·f_r625 Hz故障边带间隔Δf f_r25 Hz当齿面存在局部缺陷时每次啮合都会产生一个短暂冲击在频谱上表现为以啮合频率为中心、以转频为间隔的边带族。ToqurVibratory 的时变刚度激励本身也会产生边带这与实测中齿距误差引起的边带很难区分。因此更可靠的方法是在同一图中同时画输入轴转频处的速度幅值观察它是否与啮合频率幅值一起增大。若两者同步上升问题大概率来自轴系弯曲或不对中而不是齿面。4.4 瞬态冲击的时域判据频谱分析之外时域也可以提供快速判据。把重采样后的振动速度信号做包络然后数包络脉冲之间的时间间隔。假如间隔等于输入轴转频的倒数 1/250.04 秒说明每转一圈冲击一次指向单齿故障。如果间隔是 0.02 秒则每半圈冲击一次更多指向轴弯曲产生的双侧交替载荷。这个技巧在实践中很有效也是我对 ToqurVibratory 仿真结果做故障诊断时最先看的指标。包络提取可以直接用 MATLAB 的abs(hilbert(Vq))得到解析信号幅值然后用findpeaks找峰值位置再计算峰值间距。代码只有几行但对瞬态分析帮助极大。5. 参数修正用传递误差反推齿侧间隙的实操技巧5.1 从仿真结果计算静态传递误差齿侧间隙很难从图纸直接获得却直接影响齿轮在换向时的冲击幅度。把齿轮副的角位移按节圆半径折算到啮合线方向就能得到传递误差TE r_p·q_p - r_g·q_g。ToqurVibratory 输出的是各自由度角位移因此这一步不需要额外仿真直接对输出矩阵做线性变换。5.2 参数扫描与最小误差优化实际处理时我会把齿侧间隙作为待定参数在一定范围内扫描用实测振动响应与仿真响应之差的 RMS 做评价指标。下面是扫描 50 个间隙值的示例调用前需要在 p 结构体里补充节圆半径p.R_p、p.R_g。gaps linspace(1e-5, 1e-3, 50); rms_val zeros(50, 1); for i 1:50 p.backlash gaps(i); y0 zeros(44, 1); % 每次重新初始化 [~, Y] ode45((t, y) gear22(t, y, p), tspan, y0, opts); q Y(:, 1:22); % 传递误差 TE p.R_p * q(:, p.pinion_dof) - p.R_g * q(:, p.gear_dof); rms_val(i) rms(TE - mean(TE)); end [~, idx_min] min(rms_val); best_backlash gaps(idx_min);这段代码的关键是每次循环都重置y0。如果沿用上一次仿真的末状态当作初值系统衰减时间会叠加进评价指标导致误差曲线出现假谷值。若缺乏实测数据可以把负载平稳阶段的仿真 TE 作为参考值此时优化的是模型自洽性至少能保证参数组合不产生明显失真。5.3 验证结果与实际工况的对应反推出的齿侧间隙是否合理还需要一个快速验证。把最优间隙代回模型重新仿真后观察换向时间点的齿面瞬时接触力正常应小于额定静载荷的 2 倍。如果接触力出现断续说明间隙偏大如果完全平滑则间隙可能被刚度误差掩盖。此时可以同时扫描平均啮合刚度形成二维搜索但要注意计算量会成倍增加。用交叉验证的思路确保两个参数不会在一维曲线上同时抵消。实际操作时我会把二维搜索的结果与齿面印痕对比如果仿真接触力分布与印痕位置吻合参数才算真正可用。本文还有配套的精品资源点击获取 22自由度齿轮动力学仿真：基于MATLAB ode45的振动分析与故障诊断
简介面向机械工程与车辆工程领域的齿轮动力学学习者与研究者资源围绕22自由度齿轮传动系统的建模与数值仿真展开重点解决由齿侧间隙、时变啮合刚度等因素引发的非线性振动响应计算问题。压缩包为rar格式大小仅3KB内含1个ToqurVibratory.m脚本文件。脚本基于MATLAB平台利用ode45龙格-库塔积分器对系统状态方程进行求解涵盖了状态变量设定、质量/刚度/阻尼矩阵组装、外载荷定义、初值条件输入以及结果后处理等关键环节并用图形方式展示角位移、角速度等动态曲线。目前已有300人学习浏览。对于正在学习机械振动、齿轮动力学或MATLAB数值分析的人员这份轻量代码提供了可直接加载运行的参考模板便于深入理解非线性微分方程的求解流程并可在其基础上修改参数、扩展自由度快速开展齿轮传动系统的动态特性研究与优化验证。1. 从啮合冲击到 22 自由度ToqurVibratory 模型的切入点一台减速机齿面出现点蚀时现场振动频谱往往要到故障中后期才显露出边带。最麻烦的是你很难区分是齿形误差还是轴弯曲在激励齿轮。要回答这个问题与其对着频谱猜不如把齿轮系统建成多自由度动力学模型用仿真把啮合冲击、时变刚度和轴系扭转放到同一套方程里看。ToqurVibratory.rar 里的 ToqurVibratory.m 正是这样的实现它用 MATLAB 的 ode45 求解 22 自由度齿轮动力学方程输出角位移、角速度等状态量。这个压缩包适合传动设计、故障诊断工程师也适合想把动力学方程真正跑起来的研究生。读懂它你就掌握了一整套从建模、求解到结果的齿轮动力学工作流。2. 22 自由度齿轮模型的坐标系、自由度分配与方程组装2.1 自由度的物理含义与分配齿轮动力学计算的第一步是决定模型中哪些部件能“动”。以多级平行轴齿轮箱为例每一根轴的转动、齿轮副啮合点的相对位移、轴承座的弹性位移甚至电机转子的扭转都可能对动态响应产生影响。22 自由度意味着系统已经做成一个集中质量网络每个广义坐标对应一个弹性势能或动能项而不是把齿轮当成刚体简单折算。在 ToqurVibratory 这类模型里自由度通常按部件分组。表 2-1 是一种常见分配方式具体到某个压缩包会因齿轮级数和轴段划分而异但整体思路一致。表 2-1 22 个自由度的分组示例自由度区间物理量对应部件12角位移、横向位移电机转子 / 输入小齿轮38角位移、水平和垂直位移中间轴大齿轮、小齿轮916角位移、水平/垂直/轴向位移输出轴齿轮副与轴承座1722齿面法向相对位移、传递误差协调量多对啮合齿轮副的啮合点这里需要说明12 只是一个区间示例实际如果输入轴是一根刚性轴可能只有一个扭转自由度。但 22 自由度模型的每一列位移都会进入质量矩阵和刚度矩阵直接决定求解规模。自由度的取舍原则很简单能引起齿面载荷波动的弹性变形必须保留比如轴的扭转变形和齿轮副的啮合变形对结果影响很小的刚体平动可以合并或固定。2.2 二阶微分方程与状态变换齿轮系统的运动方程通常写成M q C q K(t) q T(t)式中 q 是 22 维广义位移向量q 和 q 是对应速度与加速度。M 是质量矩阵C 是阻尼矩阵K(t) 是时变刚度矩阵T(t) 是外载荷。K(t) 随时间变化是因为啮合点位置随齿轮转角移动齿对啮合刚度的变化是齿轮动力学区别于普通转子动力学的重要特征。ode45 要求的是常微分方程组的初值问题所以必须把二阶方程转换成一阶状态方程。常见做法是令 y [q; q]得到 44 维状态向量。下面代码展示如何在 MATLAB 函数里定义这个过程。function dydt gear22(t, y, p) % 22自由度齿轮系统的一阶状态方程 % y(1:22) 为广义位移, y(23:44) 为广义速度 q y(1:22); v y(23:44); % 时变啮合刚度: 平均刚度叠加一阶余弦波动 % p.z 齿数, q(1) 为主动轮角位移 k_m p.k_mean p.k_amp * cos(p.z * q(1) p.phi0); % 刚度矩阵中啮合位置耦合项 K p.K0; K(p.pinion_dof, p.gear_dof) K(p.pinion_dof, p.gear_dof) - k_m; K(p.gear_dof, p.pinion_dof) K(p.gear_dof, p.pinion_dof) - k_m; K(p.pinion_dof, p.pinion_dof) K(p.pinion_dof, p.pinion_dof) k_m; K(p.gear_dof, p.gear_dof) K(p.gear_dof, p.gear_dof) k_m; % 外力: 输入扭矩, 负载扭矩, 啮合误差激励 T zeros(22, 1); T(p.pinion_dof) p.T_in; T(p.gear_dof) -p.T_load; % M q T - C v - K q accel p.M \ (T - p.C * v - K * q); dydt [v; accel]; end这段代码的逻辑是先取出位移和速度再根据当前角度更新啮合刚度。刚度矩阵中同一组啮合自由度会出现 k_m 和 -k_m本质上是把两齿轮通过弹簧连接起来。p.pinion_dof和p.gear_dof是预先存好的自由度编号比如 2 和 3这样写比直接填数字更不容易出错。这里的参数p.k_mean是平均啮合刚度单位 N/mp.k_amp是刚度波动幅值通常取平均值的 10%30%。p.phi0是初始相位用来对齐齿距误差。注意M \是 MATLAB 左除它对稀疏矩阵尤其高效如果 M 是常矩阵可以在主脚本里提前做p.invM inv(M)然后用accel p.invM * (...)但要注意逆矩阵会丢失稀疏性22 自由度规模无所谓更大模型建议保留左除。2.3 参数表与激励设置要真正让模型算得动手里的几何参数必须换算成 M、C、K。表 2-2 给出一组常见示例值便于调试时对照。表 2-2 齿轮动力学计算常用参数示例参数符号示例值说明齿数z25齿轮副主从动轮独立给平均啮合刚度k_mean8e7 N/m按 ISO 6336 查表或有限元计算刚度波动幅值k_amp1.6e7 N/m一般是平均刚度的 20%阻尼比ξ0.03钢-钢齿轮副的经验值输入扭矩T_in50 N·m电机额定输出输入转速n_in1500 rpm决定啮合周期把表里的刚度除以质量矩阵对应元素就能估算系统的固有频率。例如某个齿轮副等效质量 0.05 kg啮合刚度 8e7 N/m单自由度固有频率约为 6360 rad/s约 1012 Hz。这个频率远高于转频但低于采样率所以后处理时要用足够小的最大步长才能捕获。外载荷不一定是常数。实际工况里负载突变会造成扭转冲击ToqurVibratory 这类模型在T(22) -p.T_load的位置可以改为时间分段函数。比如模拟加载瞬间可以用p.T_load * (t 0.1)让负载在 0.1 秒后加入观察瞬态冲击响应。3. ode45 求解 22 自由度系统的步长控制与刚性问题3.1 ode45 的适用边界与方法特点ode45 是 MATLAB 内置的自适应龙格-库塔法它同时使用四阶和五阶公式通过两者差估算局部误差自动调整步长。对齿轮动力学模型啮合刚度随时间连续变化不存在剧烈的尺度分裂所以 ode45 往往比隐式方法更快。但 22 自由度系统的固有频率分散如果刚度阵里同时有 1e4 和 1e8 N/m 的数量级局部误差控制会迫使步长变小仿真时间猛增。这时首先要检查的是单位而不是怀疑求解器。常用做法是先把齿轮副的啮合频率算出来用1/(20*fm)作为最大步长上界。这样不管误差控制怎么变化高频成分至少每个周期采 20 个点。为了保证 FFT 时频率分辨率够用还要保证仿真时间至少 200 个啮合周期。3.2 ode45 的标准调用方式在 MATLAB 主脚本中调用 gear22 的方式如下。这段代码可以直接替换 ToqurVibratory.m 里的求解部分前提是你有对应的 p 参数结构体。% ---- 主求解脚本 ---- p.M 0.05 * eye(22); % 演示用对角质量阵, 实际按转动惯量计算 p.C 0.03 * eye(22); % 比例阻尼, 单位 N·s/m p.K0 5e6 * eye(22); % 基础刚度阵, 单位 N/m p.k_mean 8e7; % 啮合平均刚度 p.k_amp 1.6e7; % 啮合刚度波动 p.z 25; % 主动轮齿数 p.phi0 0; % 初始相位 p.pinion_dof 2; % 主动轮自由度编号 p.gear_dof 3; % 从动轮自由度编号 p.T_in 50; % 输入扭矩 N·m p.T_load 50; % 负载扭矩 N·m y0 zeros(44, 1); % 初始静止 tspan [0 1]; % 仿真 1 秒 opts odeset(RelTol, 1e-7, AbsTol, 1e-9, MaxStep, 5e-5); [t, Y] ode45((t, y) gear22(t, y, p), tspan, y0, opts);这里的关键是AbsTol设置成 1e-9因为齿轮位移量级很小尤其是齿面法向变形可能只有微米级如果只调RelTol绝对误差过大会把微变形直接抹平。MaxStep5e-5在 1500 rpm、齿数 25 的工况下对应啮合频率 625 Hz步长 5e-5 秒时每个啮合周期约 32 步足够看清主要振动。3.3 求解精度与效率的平衡技巧表 3-1 给出 ode45 常用参数推荐这个表也适用于其他刚度和质量量级的齿轮模型。表 3-1 ode45 求解参数推荐参数推荐范围调整依据RelTol1e-6 ~ 1e-8需要捕捉边带时收紧AbsTol1e-8 ~ 1e-10位移量级为微米时取 1e-10MaxStep1/(20fm) ~ 1/(50fm)啮合频率高时取小值InitialStep自动大部分情况不需要手动工程上我的习惯是先用RelTol1e-6跑通再收紧到 1e-8比较两次结果的幅值差。如果幅值差超过 5%说明模型可能进入了刚性问题区域。判断方法很简单把opts.Stats打开看 ode45 是否频繁拒绝步长。如果 rejected 数量达到 accepted 数量的一半以上就该考虑换 ode15s。3.4 刚性问题排查与求解器切换齿轮动力学中的刚性问题通常来自轴承油膜刚度或齿面接触刚度的数值量级差异。比如齿面刚度 1e8 N/m而扭振轴段刚度 1e5 N·m/rad这会让雅可比矩阵特征值相差 3 个数量级。ode45 为了稳定会不断减小步长最终表现为卡在某个时间点不动。遇到这种情况先做三件事。第一把质量矩阵和刚度矩阵各元素的单位统一确认没有 N 和 N·m 混用。第二把阻尼比调大到 0.05 再试增加阻尼可以平滑瞬态但不解决本质问题。第三直接切到ode15s只需要把ode45换成ode15s其余代码不动[t, Y] ode15s((t, y) gear22(t, y, p), tspan, y0, opts);ode15s 是隐式变阶求解器适合刚性系统但它每个步长要解线性方程组22 自由度规模下代价很低。切过去后要重新验证时间序列是否与 ode45 在小步长下收敛到同一曲线。通常两者在驱动转矩恒定且不存在突变时差异很小如果波动很大往往不是求解器的问题而是刚度矩阵里出现了额外自由度冲突比如两个齿轮自由度重复表达了同一个位移。4. 从时间序列到频谱边带ToqurVibratory 结果的后处理流程4.1 变步长数据的等间隔采样ode45 返回的 t 和 Y 是非等间隔的。要计算 FFT必须先把信号重采样到等间隔时间轴。这里我一般选fs 10000Hz能覆盖啮合频率及其高次谐波。用interp1插值时要避免spline产生过冲pchip更适合振动信号。% 取第 8 个自由度的横向振动速度 fs 10000; t_end t(end); tq 0:1/fs:t_end; Vq interp1(t, Y(:, 822), tq, pchip); % 注意索引, 位移22 才是速度 % 去除直流分量, 便于观察振动幅值 Vq Vq - mean(Vq); figure; plot(tq, Vq * 1000); xlabel(时间 (s)); ylabel(振动速度 (mm/s));这里取的是Y(:, 822)因为前 22 列是位移后 22 列是速度。Vq单位是 m/s乘以 1000 得到 mm/s方便与现场测点数据比较。如果看的是角速度单位是 rad/s后处理时不需要换算直接做频谱即可。4.2 频率谱计算与啮合频率定位等间隔采样后用 fft 计算单边幅值谱。下面这段代码会去掉镜像部分并只显示前 3000 Hz 的频率成分一般足够覆盖齿轮故障相关频带。L length(Vq); Yf fft(Vq); P2 abs(Yf / L); P1 P2(1:L/21); P1(2:end-1) 2 * P1(2:end-1); f fs * (0:(L/2)) / L; % 只显示 3000 Hz 以内 idx f 3000; figure; plot(f(idx), P1(idx)); xlabel(频率 (Hz)); ylabel(幅值);注意P1(2:end-1)2*P1(...)是单边谱的标准做法直流分量不需要加倍。仿真数据没有传感器噪声所以频谱会比实测干净很多偶尔会出现非常尖锐的谱线这是正常现象。实际分析时要把频率分辨率dffs/L计算出来如果df大于 5 Hz可能需要延长仿真时间或补零来提高谱线密度。4.3 齿轮振动特征频率与边带识别在齿轮动力学仿真中最常看的三个频率成分是轴转频、啮合频率和边带。表 4-1 列出它们的计算方式。表 4-1 齿轮动力学特征频率计算特征公式示例值(1500rpm, z25)输入轴转频f_r n/6025 Hz啮合频率f_m z·f_r625 Hz故障边带间隔Δf f_r25 Hz当齿面存在局部缺陷时每次啮合都会产生一个短暂冲击在频谱上表现为以啮合频率为中心、以转频为间隔的边带族。ToqurVibratory 的时变刚度激励本身也会产生边带这与实测中齿距误差引起的边带很难区分。因此更可靠的方法是在同一图中同时画输入轴转频处的速度幅值观察它是否与啮合频率幅值一起增大。若两者同步上升问题大概率来自轴系弯曲或不对中而不是齿面。4.4 瞬态冲击的时域判据频谱分析之外时域也可以提供快速判据。把重采样后的振动速度信号做包络然后数包络脉冲之间的时间间隔。假如间隔等于输入轴转频的倒数 1/250.04 秒说明每转一圈冲击一次指向单齿故障。如果间隔是 0.02 秒则每半圈冲击一次更多指向轴弯曲产生的双侧交替载荷。这个技巧在实践中很有效也是我对 ToqurVibratory 仿真结果做故障诊断时最先看的指标。包络提取可以直接用 MATLAB 的abs(hilbert(Vq))得到解析信号幅值然后用findpeaks找峰值位置再计算峰值间距。代码只有几行但对瞬态分析帮助极大。5. 参数修正用传递误差反推齿侧间隙的实操技巧5.1 从仿真结果计算静态传递误差齿侧间隙很难从图纸直接获得却直接影响齿轮在换向时的冲击幅度。把齿轮副的角位移按节圆半径折算到啮合线方向就能得到传递误差TE r_p·q_p - r_g·q_g。ToqurVibratory 输出的是各自由度角位移因此这一步不需要额外仿真直接对输出矩阵做线性变换。5.2 参数扫描与最小误差优化实际处理时我会把齿侧间隙作为待定参数在一定范围内扫描用实测振动响应与仿真响应之差的 RMS 做评价指标。下面是扫描 50 个间隙值的示例调用前需要在 p 结构体里补充节圆半径p.R_p、p.R_g。gaps linspace(1e-5, 1e-3, 50); rms_val zeros(50, 1); for i 1:50 p.backlash gaps(i); y0 zeros(44, 1); % 每次重新初始化 [~, Y] ode45((t, y) gear22(t, y, p), tspan, y0, opts); q Y(:, 1:22); % 传递误差 TE p.R_p * q(:, p.pinion_dof) - p.R_g * q(:, p.gear_dof); rms_val(i) rms(TE - mean(TE)); end [~, idx_min] min(rms_val); best_backlash gaps(idx_min);这段代码的关键是每次循环都重置y0。如果沿用上一次仿真的末状态当作初值系统衰减时间会叠加进评价指标导致误差曲线出现假谷值。若缺乏实测数据可以把负载平稳阶段的仿真 TE 作为参考值此时优化的是模型自洽性至少能保证参数组合不产生明显失真。5.3 验证结果与实际工况的对应反推出的齿侧间隙是否合理还需要一个快速验证。把最优间隙代回模型重新仿真后观察换向时间点的齿面瞬时接触力正常应小于额定静载荷的 2 倍。如果接触力出现断续说明间隙偏大如果完全平滑则间隙可能被刚度误差掩盖。此时可以同时扫描平均啮合刚度形成二维搜索但要注意计算量会成倍增加。用交叉验证的思路确保两个参数不会在一维曲线上同时抵消。实际操作时我会把二维搜索的结果与齿面印痕对比如果仿真接触力分布与印痕位置吻合参数才算真正可用。本文还有配套的精品资源点击获取