直齿轮时变啮合刚度是齿轮动力学里最基础的激励源用势能法在Matlab里把它算出来顺带把齿间摩擦力也装进模型这件事网上教程不少但能一次性跑通并且拿到合理曲线的很少。我自己从推导公式到写代码、调参数折腾了一周多踩了不少坑最后把整套程序整理成了可复用的求解模型。这篇就从头到尾拆一遍为什么刚度是时变的、势能法每一步怎么推、摩擦力怎么加进模型里、Matlab代码怎么组织以及哪些地方最容易出错。如果你正在做齿轮传动系统动力学仿真或者写论文需要时变啮合刚度曲线甚至只是想搞明白能量法这东西到底怎么落地这篇文章应该能帮你省掉几天的搜索和试错时间。1. 为什么直齿轮的啮合刚度是时变的很多人第一次接触时变啮合刚度这个词会有点懵齿轮都是刚性的刚度怎么会随时间变其实这里的刚度指的是啮合齿对抵抗法向变形的能力它随啮合点位置和参与啮合的齿数交替而周期变化。1.1 单双齿交替啮合是根本原因直齿轮传动时重合度通常是个 1 到 2 之间的数比如常见的 1.4、1.6。这意味着在一个啮合周期里大部分时间是两对齿同时参与承载双齿区但有那么一小段时间只有一对齿在扛活单齿区。双齿区刚度是两个齿对刚度的并联总量大单齿区只剩一个齿对总量明显掉下来。所以刚度曲线天然就是一个周期性的锯齿波或方波圆弧过渡形状频率就是啮合频率。这个交替过程就是齿轮系统最主要的内部激励源。无论做振动响应、故障诊断还是噪声分析都得先把这条刚度曲线算对。用有限元可以算但是每换一个啮合位置就要重新建模或者重新剖分自由度大、耗时长势能法的优势在于用解析积分直接算几个刚度分量几十毫秒就能扫完一个完整啮合周期做参数扫掠、多目标优化都很合适。1.2 接触点沿齿面移动带来的几何效应除了单双齿交替还有一个容易被忽略的变化啮合点沿齿廓从齿根往齿顶移动。轮齿可以看作是固接在轮体上的变截面悬臂梁载荷作用点越靠近齿顶力臂越长弯曲变形越大等效刚度越小载荷作用点靠近齿根力臂短刚度就大。因此即便始终是单齿啮合刚度也会随着转角变化。这两个因素叠加起来最终得到的刚度曲线才有明显的碗形特征——每个单齿区像一个下凹的谷底谷底内部的细微变化来自接触点位置移动谷与谷之间的平台则是双齿区的高刚度段。理解了这一点写程序和后处理的时候心里就有数了。2. 势能法的理论骨架四个刚度分量怎么推势能法的思想说起来很简单把啮合力做的功看成几部分变形势能的叠加每一部分对应一个等效刚度最后串联得到单齿对的总刚度。这里的势能包括弯曲势能、剪切势能、轴向压缩势能和赫兹接触势能四部分。2.1 把轮齿看成变截面悬臂梁弯曲刚度以直齿轮的单个轮齿为对象把它简化为固支在齿根圆处的变截面悬臂梁。齿面法向啮合力F在啮合点作用于齿面分解为水平分量和轴向分量。水平分量使齿发生弯曲和剪切变形轴向分量使齿发生压缩变形。弯曲刚度k_b由弯曲应变能导出。设啮合点到齿根固定端的距离为d在离齿根距离为x的截面上截面面积矩决定该处的抗弯惯性矩I_x那么总的弯曲柔度可以写成积分形式1/k_b ∫₀ᵈ (F_b·(d-x))² / (E·I_x) dx这里F_b是水平分力E是弹性模量。注意I_x不是常数因为齿厚随半径变化。直齿轮齿廓是渐开线把齿厚S_x随半径r的关系写出来惯性矩I_x S_x³·L/12L是齿宽。积分从齿根x0一直积到啮合点xd每一步的I_x都不一样所以必须数值积分。2.2 剪切刚度与轴向压缩刚度剪切变形用剪切势能等效公式是1/k_s ∫₀ᵈ (1.2·F_b²) / (G·A_x) dx其中G E/[2(1ν)]是剪切模量A_x S_x·L是截面面积系数1.2是矩形截面的剪切修正系数。单看这一项的柔度占比会发现它比弯曲项小不少实测大概只有弯曲项的1/5到1/10但如果不加算出来的刚度曲线会整体偏高几个百分点在需要精确对比的实验场合不能省。轴向压缩刚度更简单只有一根杆的压缩变形1/k_a ∫₀ᵈ F_a² / (E·A_x) dxF_a是轴向分力。因为F_a本身就比较小这项柔度占比很低通常占不到总柔度的百分之几。不过既然写程序了几行积分的事一起算上既严谨又不费事。2.3 赫兹接触刚度与整体串联公式两齿面接触处的局部弹性压扁用赫兹接触理论处理接触柔度近似为常数1/k_h 4(1-ν²) / (π·E·L)注意这里的E按两轮材料参数折算成等效弹性模量。钢对钢啮合时千万别直接把206GPa代进去应当用E* 1 / [(1-ν₁²)/E₁ (1-ν₂²)/E₂]钢对钢大约113GPa左右。这也是网上很多程序算出来刚度普遍偏高一个档次的原因之一。四个分量各自代表一种串联变形路径所以单对齿的总柔度是四者之和1/k_t 1/k_b 1/k_s 1/k_a 1/k_h求逆后得到单对齿的啮合刚度k_t。如果当前处于双齿啮合区还要把正在啮合的两对齿刚度并联即k_total k_t1 k_t2。这里特别注意并联是在刚度层面相加不是柔度相加。刚接触这一块的人经常在这里犯迷糊把倒数相加得到完全错误的数量级。3. 几何前置量的精确计算啮合线与重合度程序要扫出整条时变刚度曲线最先要解决的不是刚度本身而是一堆几何量基圆半径、齿顶圆半径、标准压力角、任意半径处的齿厚、啮合线长度、重合度以及每个啮合位置对应到齿廓上的哪个点。这些量错了后面的积分全是白算。3.1 基本几何参数与渐开线齿厚公式直齿轮标准参数模数m、齿数z、压力角α、齿顶高系数ha*、顶隙系数c*。基圆半径r_b r·cosαr mz/2是分度圆半径齿顶圆半径r_a r ha*·m齿根圆半径r_f r - (ha*c*)·m。任意半径r_x处的齿厚公式用渐开线函数表示S_x 2·r_x·(π/(2z) invα - invα_x)其中invα tanα - αinvα_x tanα_x - α_x而α_x arccos(r_b/r_x)是该半径处的渐开线压力角。这个公式我建议写到程序里之前先手工验算一下在分度圆r_xr时S_x应当等于πm/2代入公式确实如此在齿顶圆处S_x小于分度圆齿厚在齿根圆处S_x大于分度圆齿厚这是自检的方向。我代码里专门写了一个assert防止inv项的符号写反。3.2 重合度计算与单双齿区域划分有了基圆和齿顶圆参数重合度按标准公式算ε [√(r_a1²-r_b1²) √(r_a2²-r_b2²) - C·sinα] / (π·m·cosα)其中C m(z1z2)/2是中心距。算出来如ε1.6说明在一个啮合周期内有60%的时间是双齿啮合40%的时间是单齿啮合。用程序模拟时我习惯把主动轮转角作为自变量把啮合线长度分成若干个等间距的离散点每个点对应一个啮合位置。区域判断的逻辑是以某对齿刚进入啮合为起点啮合线长度为g从起点到g-length对应双齿区然后进入单齿区当下一对齿开始啮合时又回到双齿区。更简单的处理办法是把整个啮合线长度等分成n个离散点预计算每个点的位置索引凡是落在双齿区窗口内的点同时计算两对齿的啮合参数。判断双齿区其实就一句话当前离散点到啮合起点的距离是否小于(p_b·(ε-1))这里的p_b πm·cosα是基圆齿距。这个条件我用在一次跑出曲线之后才发现之前用if循环判断每次都对不齐换成这个直接的关系式就干净多了。3.3 啮合点到齿根的距离d怎么取每对齿在某个啮合点的等效悬臂梁长d严格说应该从啮合点沿齿高方向投影到齿根固定端。工程上常用的简化处理是d取啮合点半径r_c与齿根圆半径r_f的差值即d r_c - r_f。当啮合点从齿顶移动到齿根时d从最大值逐渐减小到0附近物理意义直观程序也好写。如果需要更严格一点可以把齿根过渡圆弧引起的截面变化也建模进去但对大多数齿轮动力学分析来说d r_c - r_f的精度已经足够了。你要是在跟有限元结果做对比发现刚度趋势对但绝对数值偏了几个百分点通常问题出在齿根过渡部分和轮体柔性被忽略上而不是d的取法那一步的误差只是局部性的。4. 齿间摩擦力进入势能模型的关键修正标题里特意强调齿间摩擦力也有考虑进去这说明它不是常规势能法程序的默认配置。事实上绝大多数文献里的时变啮合刚度计算都忽略摩擦因为摩擦力相对法向力小一个量级对弯曲变形的贡献看起来不显著。但如果研究的是齿轮啸叫、摩擦激励或齿面微观弹流润滑的耦合效应摩擦力对刚度曲线的非对称性影响就不能再忽略了。4.1 摩擦力改变的不只是力平衡还有弯矩摩擦力沿齿面切向作用大小为F_f μ·F_Nμ是齿面摩擦系数F_N是法向载荷。它有两个效果一是直接叠加一个切向载荷分量改变剪切项和轴向项里的受力二是对轮齿固定端产生一个附加弯矩和法向力产生的弯矩叠加进而改变弯曲刚度。也就是说摩擦力修正不只是往公式里加一项μF_N那么简单它会影响所有涉及载荷分量的积分项。以单齿受力为例设啮合点压力角为α_x法向力F分解为水平分量F_b F·cosα_x和轴向分量F_a F·sinα_x。加入摩擦后水平方向的等效弯曲载荷变成F_b F·(cosα_x ± μ·sinα_x)轴向载荷变成F_a F·(sinα_x ± μ·cosα_x)符号取决于主从动轮判定和当前处于啮入段还是啮出段。核心规律是摩擦力方向始终与齿面相对滑动方向相反而在啮合节点附近相对滑动速度会反向所以摩擦力的符号在节点两侧发生突变。4.2 摩擦系数与算法中的符号处理实际钢齿轮油润滑时的齿面摩擦系数通常在0.03到0.08之间特殊润滑工况也可能到0.1以上。在程序里把μ设成一个可修改的输入参数即可方便后续做摩擦系数的参数化研究。符号问题的处理我推荐用一个方向变量flag主动轮啮入段取啮出段取-然后直接乘到μ项上。程序里写一个嵌套函数根据啮合点相对节点的位置返回符号这样既不会漏也方便以后改成随滑动速度变化的瞬时摩擦系数。如果暂时只想先加一个固定摩擦系数验证模型那就在输入文件里把μ设成0.05跑出来看看刚度的改变趋势——实测结果是有摩擦和无摩擦的曲线在大趋势上一致但在靠近啮入和啮出两端会出现明显的不对称性这也是判断摩擦修改是否生效的快速检查点。5. Matlab程序的模块结构与关键代码整套程序我分成了三个文件主脚本负责参数定义、角度扫描循环和结果绘图一个子函数计算单对齿在指定啮合半径处的四项刚度另外一个子函数输入几何参数预计算啮合线长度、重合度和单双齿切换点。这样模块化之后换齿轮副参数只需要改主脚本开头的几行常量画图代码完全不用动。5.1 主循环与啮合位置扫描以主动轮转角为自变量假设一个完整的啮合周期对应一个基圆齿距的啮合线长度将这个周期等分成N份N取100到300都可。我用N200曲线已经足够光滑。每个转角下当前啮合点沿啮合线移动径向半径r_c由几何关系推出然后调用单齿刚度子函数。如果当前位置落在双齿区就额外计算下一对齿在对应啮合位置的刚度两个刚度相加。核心主循环示意如下% 主参数 m 3; z1 20; z2 40; alpha 20*pi/180; L 20; E1 206e3; nu 0.3; mu 0.05; % 单位N/mm % 几何 r1 m*z1/2; r2 m*z2/2; rb1 r1*cos(alpha); rb2 r2*cos(alpha); ra1 r1 m; ra2 r2 m; rf1 r1 - 1.25*m; rf2 r2 - 1.25*m; % 重合度等 pb pi*m*cos(alpha); g sqrt(ra1^2-rb1^2) sqrt(ra2^2-rb2^2) - (r1r2)*sin(alpha); eps_alpha g/pb; % 扫描 N 200; k_total zeros(N,1); theta_array linspace(0, 1, N); % 归一化啮合周期 for i 1:N s theta_array(i)*g; % 当前点离啮合起点的弧长 if s pb*(eps_alpha-1) k_total(i) pair_stiffness(s) pair_stiffness(s pb); else k_total(i) pair_stiffness(s); end end这里pair_stiffness的输入是啮合点沿啮合线的位置程序内部再把弧长换算成该齿对的实际啮合半径再代入第2节那套积分公式。把重合度写成εp_b·(ε-1)的条件判断是我反复调试后最简洁的一种写法。5.2 单齿刚度子函数与数值积分单齿刚度子函数的核心是匿名函数形式的柔度积分。因为I_x和A_x都随半径连续变化而每个离散啮合位置都需要重新计算一条被积函数曲线如果每次循环都去构造一堆匿名函数速度会慢不少。我的做法是把x方向积分变量换成r方向用integral函数直接积分。示例代码片段function kt pair_stiffness(s, geo, mat) % 根据s计算当前轮齿啮合半径r_c % 计算齿根到啮合点的距离d和齿厚函数S_r fb F*(cos(alpha_c) sg*mu*sin(alpha_c)); fa F*(sin(alpha_c) - sg*mu*cos(alpha_c)); inv_bend integral((r) (fb*(d-(r-rf))).^2./(E*S_r(r).^3*L/12), rf, rc, ArrayValued, true); inv_shear integral((r) 1.2*fb^2./(G*S_r(r)*L), rf, rc); inv_axial integral((r) fa^2./(E*S_r(r)*L), rf, rc); inv_h 4*(1-nu^2)/(pi*E_star*L); kt 1/(inv_bend inv_shear inv_axial inv_h); end注意单位体系要统一我全程用N和mm弹性模量输入206e3 N/mm²这样刚度直接就是N/mm画图时通常再除个1000变成N/μm曲线数值更符合工程习惯。如果用国际单位Pa和m也不是不行但齿厚和半径都要换算成m多一层步骤就多一点出错机会。5.3 曲线输出与自检画出刚度曲线后第一件事不是看绝对值而是看形状双齿区平台高度是否基本稳定、单齿区有没有明显的对称凹陷、齿数比引起的周期特征是否合理。我建议顺手画两个图一个是刚度的时变曲线另一个是在同一张图上叠加无摩擦模型的曲线做对比这样摩擦的影响一眼就能看出来。6. 算例结果与实际调试中踩过的坑我用一组标准参数做了验证m3z120z240压力角20°齿宽20mm钢齿轮材料参数。重合度算出来约1.62单齿周期约占全周期的38%。刚度曲线呈现典型的双齿高台-单齿低谷交替形态双齿区大致在240N/μm附近单齿区掉到120N/μm左右考虑摩擦后曲线在啮入端和啮出端出现轻微不对称趋势和文献吻合。这类数量级在工程上是合理的。你可以拿这个大致范围检验自己的程序如果算出来单齿刚度大得离谱比如上千N/μm多半是等效弹性模量用错了、或者齿厚公式中inv项符号反了如果刚度曲线完全没有高低交替那大概率是单双齿判断逻辑写错了。6.1 单位不统一导致的离奇数量级这是我第一次跑出10⁹ N/mm这种荒谬结果的首因。模数用了mm弹性模量却用了PaN/m²最后所有长度项混在一起数量级差了三个度。建议在程序前面统一声明所有长度单位mm所有力单位N弹性模量单位N/mm²并写个简短注释。这类错误不调试一次根本想不到。6.2 齿厚公式里inv渐开线函数的符号渐开线函数invα tanα - α在标准计算中任意半径处齿厚要比较invα和invα_x的大小。因为齿顶处压力角大invα_x也大所以分度圆齿厚公式里是(π/(2z) invα - invα_x)。这个式子里invα_x前的负号一旦写错算出来的S_x从齿根到齿顶不降反升直接导致弯曲刚度积分异常。我的自检方法很简单把分度圆半径代进去看S_x是不是恰好等于πm/2再在齿顶处验证S_x确实变小。6.3 双齿区总刚度忘掉并联计算双齿区时有两对齿同时承载总刚度应该是两对齿各自刚度的并联k_total k1 k2。这里的k1和k2都是单对齿总刚度(已经串联了四部分柔度)。我第一次实现时下意识把两个单齿柔度相加再取倒数等于把两个弹簧串联了得到的双齿区刚度反而比单齿区还小曲线形状完全反了。这个问题只要记住多对齿同时接触是分配载荷并联加大刚度就永远不会再犯。6.4 摩擦方向突变引起的曲线毛刺考虑摩擦力之后啮合节点两侧摩擦力方向反转。用固定的号跑完整个周期时会在节点对应的位置看到一个小台阶或毛刺这不是积分误差而是物理方向突变。处理方式有两个一是明确告诉自己在节点处做符号切换二是在后处理画图时标注出节点位置不要误以为是算法不稳定而去平滑它。如果后续要分析二阶量这个突变反而是理解摩擦激励频率成分的重要信息。6.5 积分精度和被积函数的高频振动齿根附近齿厚小惯性矩I_x是三次方项被积函数在靠近齿根处变化剧烈。我用integral的默认容差误差比较大特别是齿轮模数小、齿数少时积分结果会有毫米级的微小波动。建议给integral显式设置RelTol为1e-6代价是单次扫描多花几十毫秒但对整条曲线的光滑度和可信度帮助很大。反正N200的扫描总共也只有几百次积分几十毫秒的代价完全可以接受。整套程序目前对标准直齿轮、非变位齿轮非常稳定如果后面想扩展到斜齿轮或变位齿轮还需要在几何部分增加螺旋角投影和变位系数的齿厚修正积分框架本身不用大改。我个人在跑通这个模型之后最大的体会是势能法是一个非常透明的方法——每一步积分都有明确的物理对应哪里出了问题都能回溯到具体参数。这也正是我后来在做参数敏感性分析时宁愿多花点时间用势能法、也不太愿意直接套有限元黑盒的原因。
