简介齿轮传动因承载大、效率高在机械系统中广泛应用而其动态性能受时变啮合刚度显著影响。对于斜齿轮啮合点沿齿面移动导致刚度周期性变化精确计算该参数是动力学分析与设计优化的关键。该MATLAB源程序围绕斜齿齿轮时变啮合刚度计算展开通过定义模数、压力角、螺旋角、齿宽等基本几何参数结合弹性变形模型与接触状态分析完成局部刚度及整体啮合刚度的求解并提供可视化展示与频域处理辅助函数。资源包共2个文件均为.m源码文件大小仅1KB包含主程序与数据处理脚本轻量易读适合机械、航空航天及汽车领域的工程师和研究人员参考。已有950人学习下载可用于理解齿轮啮合刚度算法流程、快速运行代码并观察刚度波动特性为齿轮系统振动噪声预测与传动优化提供实用工具。1. 时变啮合刚度为什么是齿轮动力学仿真的第一道坎齿轮箱的振动噪声、齿面点蚀剥落、传动误差预测最后都会收敛到一个共同的激励源啮合刚度。很多人第一次拿到“calculation of gear mesh stiffness.zip”这个MATLAB包时以为只是把齿轮几何参数代进公式就能得到一条直线实际上齿轮啮合过程中单双齿交替、接触点沿齿廓滑动、轮体弹性变形同时发生刚度是一个随啮合位置周期性波动的曲线波动的幅度和相位直接决定动态响应。做齿轮箱故障诊断的工程师如果拿恒定刚度去算边频带的来源根本解释不通做传动系统设计的如果只看均值共振转速的预测会整体偏移。这篇文章把时变啮合刚度从理论模型讲到MATLAB代码实现再到结果验证和工程扩展适合要做动力学仿真、参数扫参或者基于振动信号的齿轮检测算法开发的人。2. 势能法计算齿轮时变啮合刚度的理论模型与关键参数2.1 势能法把单齿对啮合刚度拆成五个串联分量时变啮合刚度的计算思路有很多条路有限元法最接近真实但网格和接触设置代价高每改一次齿轮参数就要重新建模石川法把齿廓简化为梯形梁计算快但精度粗糙目前工程和科研里用得最多的是势能法。势能法的基本假设是一对啮合齿在载荷作用下的总弹性变形可以分解为几个独立的部分——齿的弯曲变形、剪切变形、轴向压缩变形、轮体弹性变形以及齿面接触处的Hertz局部变形。这五部分在同一载荷路径上依次发生所以是串联关系1/k_pair 1/k_b 1/k_s 1/k_a 1/k_f 1/k_h其中k_b是弯曲刚度k_s是剪切刚度k_a是轴向压缩刚度k_f是轮体变形刚度k_h是Hertz接触刚度。串联的含义很直接总柔度等于各部分柔度之和刚度取倒数。把齿看作一个变截面悬臂梁从齿根到接触点沿齿高方向积分就能得到各项刚度。势能法的优势在于每个分量都有明确的几何意义齿轮参数变化后可以立刻反映到积分结果里单次计算在毫秒级非常适合后续做参数扫参和优化。2.1.1 悬臂梁模型里弯曲、剪切与压缩的积分式对于主动轮或从动轮中的任意一个齿轮设接触点到齿根的距离为d沿齿高方向取微元dx该截面的齿厚为s(x)截面惯性矩I(x) b·s(x)³/12截面积A(x) b·s(x)。接触点压力角为α_k法向力F分解为水平分量F·cosα_k和垂直分量F·sinα_k则弯曲、剪切、轴向压缩的柔度分别为1/k_b ∫₀ᵈ [cosα_k·(d−x) sinα_k·h_c]² / (E·I(x)) dx1/k_s ∫₀ᵈ 1.2·sin²α_k / (G·A(x)) dx1/k_a ∫₀ᵈ cos²α_k / (E·A(x)) dx其中h_c是载荷偏心距工程上常取接触点的半齿厚用于近似垂直分量产生的附加弯矩1.2是矩形截面剪切系数G E / (2(1ν))。这套公式的精度取决于齿厚函数s(x)是否准确而s(x)必须由渐开线方程给出不能用梯形近似糊弄。2.2 啮合周期、基节与重合度先算对几何再谈刚度时变啮合刚度的“时变”本质来自重合度。齿轮连续啮合时有时一对齿单独承载有时前后两对齿同时承载。单齿区刚度低、变形大双齿区刚度高、载荷被分担二者交替出现就形成了刚度波动。两个关键几何量必须在计算刚度前先确定基节p_b π·m·cosα₀以及实际啮合线长度L_a √(r_a1²−r_b1²) √(r_a2²−r_b2²) − a·sinα_w。重合度ε L_a / p_b。标准直齿轮的重合度通常在1.2到1.8之间这意味着啮合过程中有长度为(ε−1)·p_b的双齿区出现在啮入端和啮出端中间长度为p_b的单齿区。如果重合度算错双齿叠加的位置就会错位整个刚度曲线的波形都会变形这是新手最容易踩的坑。2.3 单位制与坐标约定N/mm 还是 N·m/rad势能法计算时全部采用N和mm单位制模数、齿宽、半径用mm弹性模量E用MPa即N/mm²齿宽b用mm最后得到的刚度单位是N/mm。若想换算成扭转刚度N·m/rad需要乘上半径平方k_θ k·r²。很多MATLAB程序跑出来结果量级不对十有八九是E的单位用了Pa没除以10⁶或者齿宽忘了乘进去。另一个常见约定是啮合线上位置ξ的取向我一般定义ξ0为啮入点主动轮齿顶圆与啮合线交点ξL_a为啮出点这样单双齿区的判断最直观。下表汇总五个刚度分量的量级和主要影响因素方便实现后对照检查刚度分量典型量级b20mmm3主要影响因素弯曲刚度 k_b5×10⁴ ~ 2×10⁵ N/mm模数、齿数、接触点位置剪切刚度 k_s1×10⁵ ~ 5×10⁵ N/mm齿宽、齿根厚度轴向压缩 k_a10⁵ ~ 10⁶ N/mm齿宽、压力角Hertz接触 k_h10⁶级 N/mm齿宽、弹性模量轮体变形 k_f10⁵ ~ 10⁶ N/mm齿根圆角、轮缘厚度Hertz项在总柔度中占比不大但绝对不可省略轮体变形k_f常被忽略导致总刚度高估10%~20%在齿轮检测的定量诊断里这个误差足以掩盖早期裂纹的特征。3. 用MATLAB实现齿轮时变啮合刚度计算的完整流程3.1 齿轮基本参数定义与渐开线齿厚函数实现的第一步是把齿轮几何参数定义成结构体方便后续传递给子函数。我用一组标准直齿轮参数做示例模数m3主动轮齿数z124从动轮z248压力角20°齿宽20mm。MATLAB代码从参数定义开始% 齿轮基本参数单位统一为 mm 和 MPa p.m 3; % 模数 mm p.z1 24; p.z2 48; % 齿数 p.a0 20*pi/180; % 分度圆压力角 rad p.b 20; % 齿宽 mm p.x1 0; p.x2 0; % 变位系数 p.E 2.06e5; % 弹性模量 MPa p.nu 0.3; % 泊松比 p.G p.E/(2*(1p.nu)); % 切变模量 MPa p.ha_star 1; % 齿顶高系数 p.c_star 0.25; % 顶隙系数 % 几何计算 p.r1 p.m*p.z1/2; p.r2 p.m*p.z2/2; p.rb1 p.r1*cos(p.a0); p.rb2 p.r2*cos(p.a0); p.ra1 p.r1 p.ha_star*p.m; p.ra2 p.r2 p.ha_star*p.m; p.rf1 p.r1 - (p.ha_starp.c_star)*p.m; p.rf2 p.r2 - (p.ha_starp.c_star)*p.m; p.a p.r1 p.r2; % 标准安装中心距 p.aw acos((p.r1p.r2)*cos(p.a0)/p.a); % 啮合角 p.pb pi*p.m*cos(p.a0); % 基节 mm p.La sqrt(p.ra1^2-p.rb1^2) sqrt(p.ra2^2-p.rb2^2) ... - p.a*sin(p.aw); % 实际啮合线长度 mm p.eps p.La / p.pb; % 重合度这段代码把所有半径做成结构体字段后面子函数只传p一个变量就能访问全部参数。p.eps在1.2~1.8之间是正常范围如果算出来小于1说明参数定义有问题。注意p.aw这里用的是余弦公式实际标准安装时恒等于p.a0保留这一步是为了以后改中心距变位时不用回翻代码。接下来是渐开线齿厚函数。齿廓上任意半径r处的齿厚由渐开线方程决定不能直接用梯形近似function s tooth_thickness(r, z, m, a0, x, rb) % 计算半径 r 处的渐开线齿厚r 可以是向量 r0 m*z/2; % 分度圆半径 s0 m*(pi/2 2*x*tan(a0)); % 分度圆齿厚含变位修正 ak acos(rb ./ r); % 该半径处的压力角 inva0 tan(a0) - a0; % 渐开线函数值 invak tan(ak) - ak; s 2*r .* (s0/(2*r0) inva0 - invak); end这里invak是渐开线函数的缩写即tanα−αMATLAB里没有内置函数所以要手写。变位系数x已经进了分度圆齿厚s0的表达式后续做齿轮变位优化时只需要改p.x1、p.x2。3.2 悬臂梁刚度积分的MATLAB子函数有了齿厚函数就可以沿齿高把齿离散成薄片做数值积分。我把弯曲、剪切、轴向压缩三个刚度放在一个子函数里返回避免在循环中重复计算公共量。function [kb, ks, ka] gear_tooth_stiff(rf, rk, p, z, rb) % rf: 齿根圆半径, rk: 接触点半径, p: 参数结构体 n 200; % 齿高微元数 rv linspace(rf, rk, n); % 从齿根到接触点离散 dr rv(2) - rv(1); s tooth_thickness(rv, z, p.m, p.a0, p.x1, rb); I p.b * s.^3 / 12; % 截面惯性矩 A p.b * s; % 截面积 alpha_k acos(rb / rk); % 接触点压力角 x rv - rf; % 距齿根距离 hc 0.5 * s(end); % 载荷偏心距取接触点半齿厚 Mc cos(alpha_k) * (x(end)-x) sin(alpha_k) * hc; kb 1 / sum(Mc.^2 .* dr ./ (p.E .* I)); ks 1 / sum(1.2 * sin(alpha_k)^2 .* dr ./ (p.G .* A)); ka 1 / sum(cos(alpha_k)^2 .* dr ./ (p.E .* A)); end积分离散数n取200在精度和速度之间比较平衡改成500后刚度变化通常小于0.5%。如果n取得太小比如20曲线会出现明显锯齿。Mc对应弯矩臂由水平分量的力臂(x(end)-x)和垂直分量的偏心距hc组成单位是mm平方后除以惯性矩再积分量纲正好是1/(N/mm)。3.3 单双齿啮合区拼接与啮合线位置参数化单齿对的刚度算完之后还要考虑重合度带来的双齿并联。这里用啮合线坐标ξ作为自变量最为直接。ξ从0变化到p.La主动轮和从动轮的接触点半径都可以由渐开线性质直接导出% 主动轮齿顶圆到啮合线上一点的距离决定接触半径 rk1 sqrt(p.rb1^2 (sqrt(p.ra1^2-p.rb1^2) - xi).^2); rk2 sqrt(p.rb2^2 (sqrt(p.ra2^2-p.rb2^2) - (p.La-xi)).^2);当ξ0时主动轮接触半径等于ra1齿顶从动轮接触半径最小当ξp.La时正好相反。单齿对刚度函数把五个分量串联起来function k_pair pair_stiffness(rk1, rk2, p, z1, z2, rb1, rb2) [kb1, ks1, ka1] gear_tooth_stiff(p.rf1, rk1, p, z1, rb1); [kb2, ks2, ka2] gear_tooth_stiff(p.rf2, rk2, p, z2, rb2); kh pi * p.E * p.b / (4*(1-p.nu^2)); % Hertz线接触刚度 cf 2.5e-3 * cos(p.aw)^2 / (p.E * p.b); % 轮体变形简化项 inv_k 1/kb1 1/ks1 1/ka1 ... 1/kb2 1/ks2 1/ka2 ... 1/kh cf; k_pair 1 / inv_k; end轮体变形项这里先用常系数代替严格做法是按Sainsot公式根据齿根截面几何拟合工程上把这个系数做成查表函数即可。3.4 绘制时变啮合刚度曲线的完整主程序主循环遍历整个啮合线在每个位置先算当前齿对的刚度再判断是否处于双齿区是则叠加相邻齿对的贡献N 500; xi linspace(0, p.La, N); k_mesh zeros(1, N); for i 1:N % 本齿对接触位置 rk1 sqrt(p.rb1^2 (sqrt(p.ra1^2-p.rb1^2) - xi(i)).^2); rk2 sqrt(p.rb2^2 (sqrt(p.ra2^2-p.rb2^2) - (p.La-xi(i))).^2); k_cur pair_stiffness(rk1, rk2, p, p.z1, p.z2, p.rb1, p.rb2); k_extra 0; if xi(i) (p.eps-1)*p.pb % 啮入端双齿区后一对齿在前方一个基节处 xi2 xi(i) p.pb; rk1b sqrt(p.rb1^2 (sqrt(p.ra1^2-p.rb1^2) - xi2).^2); rk2b sqrt(p.rb2^2 (sqrt(p.ra2^2-p.rb2^2) - (p.La-xi2)).^2); k_extra pair_stiffness(rk1b, rk2b, p, p.z1, p.z2, p.rb1, p.rb2); elseif xi(i) p.pb % 啮出端双齿区前一对齿在后方一个基节处 xi2 xi(i) - p.pb; rk1b sqrt(p.rb1^2 (sqrt(p.ra1^2-p.rb1^2) - xi2).^2); rk2b sqrt(p.rb2^2 (sqrt(p.ra2^2-p.rb2^2) - (p.La-xi2)).^2); k_extra pair_stiffness(rk1b, rk2b, p, p.z1, p.z2, p.rb1, p.rb2); end k_mesh(i) k_cur k_extra; end % 横坐标换成主动轮转角更直观theta1 xi / rb1 (rad) theta1 xi / p.rb1 * 180/pi; plot(theta1, k_mesh/1000, LineWidth, 1.5); xlabel(主动轮啮合转角 (deg)); ylabel(啮合刚度 (N/mm × 10^3)); grid on;双齿区的叠加逻辑是这段代码最核心的部分啮入端当前齿对刚进入啮合它前面一个基节处还有一对齿处在啮出端所以要取xi p.pb啮出端则相反取xi - p.pb。两个边界条件用p.eps-1)*p.pb和p.pb判断依据是重合度定义。如果这里写反曲线会出现一侧多峰一侧凹陷的对称错误而且均值可能严重偏离。4. 齿轮时变啮合刚度计算结果的验证方法与精度校准4.1 用ISO 6336-1经验值做数量级校核代码跑通后第一件事不是看曲线漂不漂亮而是做数量级校核。ISO 6336-1标准给出钢制标准直齿轮单齿对啮合刚度的经验参考值c大约在10~20 N/(mm·μm)换算成这里的刚度单位就是乘齿宽b再乘1000也就是2×10⁵~4×10⁵ N/mm的量级对b20mm的齿轮。我一般把MATLAB算出的单齿区曲线均值落在这个区间同时看双齿区刚度是否比单齿区高30%~70%波动幅度是否在20%~50%之间。校核项判断依据合理范围重合度 εp.La / p.pb1.2 ~ 1.8单齿区均值ISO 6336-1 经验值2×10⁵ ~ 4×10⁵ N/mm双齿/单齿刚度比理论推导1.3 ~ 1.8波动幅度 (k_max−k_min)/k_mean直齿轮典型值0.2 ~ 0.5如果偏差超过这个范围优先检查单位制而不是代码逻辑。E用了Pa而不是MPa、齿宽忘乘、刚度倒数输出时忘取倒数是三个最经典的错误点。这里的轮体变形简化项对均值影响明显调校时可以先把它置零看基线再加上看变化幅度两部分分开验证。4.2 重叠度与基节计算错误的典型表现重合度算错在曲线上的表现非常典型如果p.eps算出来大于实际值双齿区会被拉长曲线中段的单齿平台变短甚至消失整体曲线看起来“太饱满”如果算小了双齿区过短曲线出现明显的尖刺。最隐蔽的错误是基节用成了齿距πm标准直齿轮的基节是πm·cosα₀在压力角20°时两者相差约6%足以让双齿区边界偏出半个齿距位置。调试时我习惯把p.La、p.pb、p.eps这三个值先打印出来和手算值对比。另一个容易被忽视的问题是齿轮变位后齿顶圆半径变化ra1 r1 m这个式子在变位齿轮上要改成r1 (ha_star x1)*m否则重合度直接算错。做齿轮检测的同行拿这套曲线做故障特征时尤其要注意变位齿轮的边界条件否则后续包络谱分析里每条边频带都会对不上。4.3 离散微元数与积分精度改哪里先看结果曲线gear_tooth_stiff里的微元数n对结果的影响不是线性的。n从20加到200刚度均值可能上升或下降2%~5%这个方向取决于齿根附近齿厚函数的变化率n从200加到1000变化通常小于0.3%。判断是否收敛的方法是扫一遍n值画出刚度均值变化曲线而不是只看单一数值。主循环里的啮合线采样数N也有讲究N太小时双齿区边界处的刚度突变会被平滑掉导致曲线过渡区失真N取300~500对直齿轮已经足够。还有一个精度坑trapz和sum在等距网格上的积分结果几乎一致但如果不小心把dr写成rv(2)-rv(1)之后又改了n的奇偶性网格间距不匹配会引入累积误差。我的习惯是让上层调用统一传入n并且用linspace自动保证等距避免手动生成网格。5. 从单工况到批量扫参齿轮时变啮合刚度的工程扩展5.1 把MATLAB计算函数封装成批量扫参接口把第3章的主循环包成一个函数输入模数、齿数、齿宽、变位系数输出一个周期的刚度曲线和重合度就可以用循环做参数扫描。我一般会配合MATLAB优化工具箱以刚度波动幅度最小为目标约束条件限死最小齿厚和重合度下限搜索齿数比和变位系数的最优组合。多目标时注意刚度均值和波动幅度往往互相矛盾增大模数提升均值但波动幅度不一定减小扫完参数后建议把帕累托前沿画出来不要只看单目标结果。5.2 齿廓修形与齿根裂纹下的刚度扰动模拟齿顶修缘在代码层面的实现方式是在齿厚函数里对靠近齿顶的部分做厚度削减或者直接对接触点有效半径范围的上下限做截断本质上模拟的是接触线长度变短。鼓形齿修形则要引入沿齿宽方向的载荷分布系数在赫兹项里乘上一个小于1的接触系数。齿根裂纹的模拟更直接把gear_tooth_stiff里齿根到裂纹尖端的弯曲惯性矩减小裂纹越深、刚度下降越明显而且裂纹只会影响单齿区对应相位这正好是齿轮检测中识别早期剥落和断齿的物理基础。5.3 把时变刚度表接入齿轮动力学模型的查表技巧计算得到的刚度曲线最终要喂给动力学模型才有价值。常见做法是把横坐标从啮合线换成主动轮转角θ1存成两张等长的表角度、刚度然后在ode45的激励函数里用interp1(theta_table, k_table, mod(theta, T))循环查表。查表时注意啮合周期T对应主动轮转过一个基节角的弧度值不要和齿距角混在一起。如果后续要做深度学习故障诊断建议把多工况下扫出的刚度曲线归一化后存成矩阵作为特征样本输入这样比直接用原始振动信号少一层噪声干扰。最后提醒一个工程细节刚度表的采样间隔至少要覆盖啮合周期内20个点否则双齿区边界的突变在数值积分里会变成高频激励导致动力学响应里出现不真实的毛刺。按这个接口导出的刚度表可以直接作为ode45里查表激励项下一步要做的是把重合度对步长的敏感性也一并标定进去。本文还有配套的精品资源点击获取