斜齿轮时变啮合刚度计算：势能法与切片法的Matlab实现

1. 项目概述为什么斜齿轮的啮合刚度不是常数又为何这么难算做传动系统动力学和故障诊断的朋友恐怕都绕不开时变啮合刚度这个东西。它直接决定齿轮副的啮合激励是振动响应、噪声辐射、齿面动载荷分配最核心的输入参数之一。直齿轮的时变啮合刚度用经典的势能法加一个二维齿廓截面就能算得比较准但一换成斜齿轮很多人第一反应就是上有限元模型一搭就是半天网格稍微细一点还会带出收敛问题。这篇文章想聊一条折中路线用势能法处理单个齿的变形用切片法处理斜齿轮特有的螺旋接触线在Matlab里把整个模型搭起来精度足够支撑大多数工程分析速度却比有限元快两个数量级。斜齿轮的时变啮合刚度为什么难算关键在啮合线的空间形态。直齿轮接触线是一条与齿宽平行的直线任意瞬时参与啮合的齿廓截面可以当作同一个二维问题来处理。斜齿轮不一样由于螺旋角的存在齿面接触线是一条斜线从啮入到啮出不是整条齿宽同时参与而是沿着齿宽方向逐渐进入、逐渐退出。这导致同一时刻不同齿宽截面的啮合位置完全不同必须三维建模才能描述。常见的做法有三种一是完整三维有限元精度高但建模成本和计算成本都很高参数扫掠时尤其痛苦二是解析公式加修正系数速度快但适用范围窄螺旋角一大偏差就很明显三是把斜齿轮沿齿宽切成若干薄片每一片按照直齿轮的势能法计算再把所有薄片的贡献叠加起来。这个方法就是切片法也是本文程序的核心思路。它把三维问题拆成一系列二维问题兼顾精度和速度适合做参数研究和优化迭代。本文适合三类读者一是搞齿轮传动设计需要快速评估不同螺旋角、齿宽和修形参数下啮合刚度变化的工程师二是做齿轮箱振动故障诊断需要把啮合刚度作为激励源输入到动力学模型中的人三是在校学生或研究人员想把势能法和切片法从公式变成可运行程序作为后续研究的基线代码。既然是谈程序实现我就把整个建模过程拆开讲从势能法每个刚度分量的物理含义到切片法如何处理斜齿轮螺旋角再到Matlab代码的模块划分最后给出验证方法和常见翻车点。按照这个顺序读下来你应该能复现一个可以放进论文或工程报告里的程序。2. 势能法的核心拆解五个变形分量的物理对应和数学表达势能法并不是什么高深的数值方法它本质上是用材料力学的应变能原理把轮齿在外载荷作用下的变形分解成几个独立的变形模式分别求柔度再串联叠加。对直齿圆柱齿轮单个齿在啮合点受法向力后变形主要来自五部分齿面弯曲、剪切、轴向压缩、齿体弹性变形含齿根圆角、轮缘变形以及两齿面接触处的局部赫兹变形。为什么一定要拆成五个分量因为轮齿的受力状态不是简单拉压或纯弯而是悬臂梁式弯曲、剪切和局部接触同时存在把它们合并成一个刚度系数精度和物理解释都会有损失。在时变啮合刚度的计算中我们还希望看到每个分量随啮合位置的变化规律方便后续做灵敏度分析和齿形修正判断。2.1 五个柔度分量的总框架把单齿看成一根变截面悬臂梁齿根固定齿顶受载。轴向坐标沿齿高方向载荷点随啮合位置移动。按照经典势能法的形式单齿的柔度可以写成[ \frac{1}{k_t} \frac{1}{k_b} \frac{1}{k_s} \frac{1}{k_a} \frac{1}{k_f} \frac{1}{k_h} ]这里 (k_b) 是弯曲刚度(k_s) 是剪切刚度(k_a) 是轴向压缩刚度(k_f) 是齿体弹性刚度(k_h) 是赫兹接触刚度。注意这些刚度不是并联而是串联因为同一个载荷作用点后面各项变形依次发生挠度叠加对应的柔度刚度倒数相加。程序实现时我习惯先算柔度再在最后取倒数。这样每一步的数值都保持在相近量级不会因为某个刚度特别大而出现数值溢出。2.2 弯曲、剪切和压缩三个积分的实现细节弯曲柔度是五个分量里最吃计算量的一项。它的表达式通常写成对齿廓坐标的积分[ \frac{1}{k_b} \int_{0}^{d} \frac{\left[1 - \frac{x}{d} \right]^2}{E I_x} , dx ]其中 (d) 是载荷作用点到齿根固定端的距离(I_x) 是当前截面的惯性矩。对齿轮来说截面不是矩形齿廓曲线由渐开线、过渡曲线和齿根圆弧组成所以 (I_x) 必须由当前截面的齿厚计算。数值上我会先把齿廓离散成几百个点然后逐步累积积分而不是用符号积分硬推。剪切柔度表达式是[ \frac{1}{k_s} \int_{0}^{d} \frac{1.2}{G A_x} , dx ]这里的 1.2 是矩形截面剪切系数。虽然齿截面不是严格矩形但工程上大多数文献仍取 1.2误差可以接受。(G) 是剪切模量由弹性模量 (E) 和泊松比 (\nu) 换算(G E / [2(1\nu)])。轴向压缩柔度则对应齿面法向力的轴向分量[ \frac{1}{k_a} \int_{0}^{d} \frac{\sin^2 \alpha_x}{E A_x} , dx ]式中 (\alpha_x) 是当前截面处啮合力方向与齿高方向的夹角。很多初学者实现时会漏掉这一项因为直齿轮中轴向分量相对小但对斜齿轮法向力在端面的投影变化比较明显漏掉后刚度曲线会偏高。2.3 齿体弹性刚度的圆角修正齿体弹性刚度 (k_f) 描述的是齿根圆角以下、轮体部分的弹性让位。这一项不能用悬臂梁模型简单描述通常采用引入无量纲系数的经验表达式。计算时由一个关键尺寸比 (r_f / h_f) 和载荷角共同决定其中 (r_f) 是齿根圆角半径(h_f) 是轮缘厚度。实际工程中不少人会直接把齿体弹性刚度忽略尤其是轮缘较厚时。这个做法在大多数实心齿轮上是可接受的但对薄轮缘齿轮齿体弹性变形能占到总变形的百分之二三十忽略会造成明显偏差。我在程序里保留这一项并预留系数调节口方便薄轮缘工况切换。2.4 赫兹接触刚度为什么它可以当常数赫兹接触刚度反映两个齿面在接触点附近的局部弹性压陷。它的核心表达式可以写成包含接触半宽和等效曲率半径的形式但对齿轮副而言两个齿面在啮合点附近的曲率半径是随啮合位置变化的严格说也应该是变量的函数。不过工程实践中有一个非常实用的近似把赫兹接触刚度近似为与材料和接触线长度相关的常数尤其在载荷不大、齿面修形不明显时这一近似足够准确。因为赫兹刚度的数量级很大通常占总柔度比例较小它的变化对总刚度曲线的趋势影响有限。在切片法里每一片有自己的接触线长度赫兹接触刚度按比例分配到各切片上这样既保住局部性又避免过度复杂化。把这五部分组合起来单个齿截面在某啮合位置的刚度就确定了。要注意这里的 (k_t) 只是一个齿截面的刚度完整的一对齿轮啮合刚度还要同时考虑主、从动轮两个齿的变形也就是两个单齿柔度相加[ \frac{1}{k_{pair}} \frac{1}{k_{t1}} \frac{1}{k_{t2}} ]对多齿对同时啮合的情况再按并联叠加总啮合刚度等于所有同时参与啮合的齿对的刚度之和。这个逻辑在直齿轮里很清晰到了斜齿轮就要跟切片法结合这也是下一节要说的重点。3. 切片法处理螺旋接触线从直齿到斜齿的关键一步切片法的基本思想非常朴素把斜齿轮沿齿宽方向切成若干片薄片每一片近似看成直齿用上一节的势能法计算该截面在当前啮合时刻的刚度然后把所有薄片的柔度按一定规则合成。问题在于这个“一定规则”不是简单的把刚度全部加起来而是要严格还原斜齿轮的啮合时序。3.1 斜齿轮啮合的时序差异从哪来直齿轮沿齿宽看所有截面的啮合状态完全同步所以相加就是几对齿的刚度直接求和。斜齿轮有一个齿宽方向的相位差由于齿面是螺旋面同一时刻不同截面上同一对齿所处的啮合位置是不同的。可以这样理解你把齿宽方向的一条线想象成拧在圆柱上的螺旋线螺旋线每走过一段轴向距离端面上对应的角度就往前推进一段。这个角度差就是相位差。这个相位差用起来有很直观的表现斜齿轮的接触线是一条从齿顶到齿根斜着跨越齿宽的对角线实际啮合时齿的一端先接触紧接着接触区沿齿宽方向逐渐扩展到另一端然后退出。对应的刚度变化就不是直齿轮那种台阶式跳变而是相对平滑的上升和下降。这也是斜齿轮传动比直齿轮更平稳的根本原因。3.2 切片几何和相位滞后的数学描述设斜齿轮的螺旋角为 (\beta)齿宽为 (b)沿齿宽方向切成 (N) 片每片厚度为[ \Delta s \frac{b}{N} ]第 (i) 片相对于第一片的轴向距离是 (i\cdot \Delta s)对应的端面角度滞后为[ \Delta \phi_i \frac{i\cdot \Delta s \cdot \tan\beta}{r_b} ]其中 (r_b) 是基圆半径。这个滞后量非常关键。计算整条接触线上的啮合状态时不能把第 0 片的啮合位置直接套到第 (i) 片上而要在滚动角上偏移 (\Delta\phi_i)。我在初版程序里就吃过这个亏把每片的接触起始位置都当成相同算出来的斜齿轮刚度曲线非常接近直齿轮波动明显偏大。后来对照文献里类似算例才发现相位滞后遗漏了。3.3 切片内直齿计算与整条接触线的合成每一片内的计算完全复用直齿齿轮的势能法流程但要传入一个局部啮合角度。局部啮合角度由全局啮合转角加上该切片自身的相位滞后得到[ \theta_i \theta \Delta\phi_i ]这个局部角度决定该切片当前是处于啮入、完整啮合还是啮出阶段也决定该片齿面上受载点的具体位置。得到每片的单齿柔度之后同一齿对内所有参与接触的切片按串联方式合成因为同一对齿的接触载荷是流过所有参与截面的[ \frac{1}{k_{pair}} \sum_{i} \frac{1}{k_{slice,i}} ]不同齿对之间则是并联因为载荷被多对齿共同分摊总刚度等于并联刚度之和。这个串联、并联关系是切片合成里最容易搞混的地方我会在下一节代码实现里反复强调。3.4 切片数量怎么定切片数 (N) 直接决定计算精度和速度。我做过一组对照实验模数 2.5、齿宽 40mm、螺旋角 16° 的斜齿轮切片厚度从 2mm 依次减到 0.2mm总刚度的变化逐渐收敛。切片厚度超过 1mm 时刚度曲线有明显的锯齿状波动小于 0.5mm 后曲线基本重合最大偏差不到百分之三。经验取值是让切片厚度小于 0.5mm或者切片数量不少于 80。如果你的齿宽很大比如超过 200mm把切片数量提高到 200 以上保证每片宽度不超过 1mm。再密就没有必要了计算时间会线性增长精度收益却很小。4. Matlab程序实现输入参数、几何离散和刚度合成主流程到这里理论部分已经具备下面进入实际可运行的实现。我在Matlab中把程序分成三个层次主程序负责扫掠啮合转角并输出刚度曲线一个函数负责计算单片的五项刚度一个函数负责几何和派生参数计算。这样的结构适合参数扫掠也方便后续替换齿轮几何。4.1 输入参数和派生几何量主程序开头是一组输入参数以法向模数、齿数、压力角、螺旋角和齿宽为主同时给出材料和切片设置%% 基本输入参数 mn 2.5; % 法向模数单位 mm z1 23; % 主动轮齿数 z2 47; % 从动轮齿数 alpha_n 20 * pi / 180; % 法向压力角单位 rad beta 16 * pi / 180; % 分度圆螺旋角单位 rad b 40; % 齿宽单位 mm ds 0.5; % 切片厚度单位 mm %% 材料参数 E 206000; % 弹性模量单位 MPa nu 0.3; % 泊松比 %% 派生参数 mt mn / cos(beta); % 端面模数 alpha_t atan(tan(alpha_n) / cos(beta)); % 端面压力角 rb mt ./ 2 .* [z1 z2] .* cos(alpha_t); % 基圆半径 ra mt ./ 2 .* [z1 z2] mn; % 齿顶圆半径 rf mt ./ 2 .* [z1 z2] - 1.25 .* mn; % 齿根圆半径近似 Nslice round(b / ds); % 切片个数这里有几个细节要说明。端面模数是法向模数除以螺旋角的余弦因为斜齿轮的几何尺寸是在端面度量的。端面压力角也同理法向压力角不能直接用于基圆半径计算。齿顶圆半径和齿根圆半径这里给了简化表达实际工程齿轮可能有变位、修顶、齿顶倒角需要按齿轮手册的公式替换。程序结构上把这些参数放在一起就是为了后续替换方便。4.2 单齿截面的五项刚度计算函数每一个切片内我需要根据局部啮合角度计算该截面对应的受载位置再对齿高方向做数值积分。函数接口可以设计成这样function k_slice calc_slice_stiffness(theta_local, rb_now, ra_now, E, nu, mn, z, beta) % theta_local: 该切片的局部啮合转角单位 rad % rb_now: 当前轮的基圆半径mm % ra_now: 当前轮的齿顶圆半径mm % 返回: 该截面的单齿等效刚度单位 N/mm % 1. 由局部转角计算载荷点位置 % 2. 沿齿高方向离散若干截面 % 3. 对弯曲、剪切、压缩、基体柔度积分 % 4. 叠加赫兹接触柔度 % 5. 取倒数得到单齿刚度 end具体积分可以用 for 循环沿齿高方向走也可以向量化。齿数少时两者差距不大齿数多且切片多时建议向量化。我用过向量化的写法把一个齿的所有截面同时算出面积、惯性矩、柔度贡献速度快不少但代码可读性稍微下降读者按自己习惯取舍。齿廓离散点是程序的隐形关键。理论渐开线齿廓到离齿根一定距离后就不再存在实际的齿廓的过渡段由齿根圆弧或摆线构成。处理办法从齿根圆到基圆之间采用过渡圆弧近似从基圆到齿顶圆用渐开线方程离散。这个处理是否合理会直接影响齿体弹性刚度的精度也影响弯曲刚度的积分下限。4.3 整个接触线上的刚度合成单片的刚度算完后接下来就是把所有切片组装。下面这段是主循环的核心框架%% 啮合转角扫描范围 theta_step 0.02; % 角度步长单位 rad theta_vec 0 : theta_step : 2*pi; Kmesh zeros(size(theta_vec)); for j 1 : length(theta_vec) theta_global theta_vec(j); % 对当前全局转角分别计算每对齿的贡献 % 先考虑单对齿遍历所有切片 inv_k_pair 0; for i 1 : Nslice s_axial (i - 0.5) * ds; phi_shift s_axial * tan(beta) / rb(1); theta_local theta_global phi_shift; % 判断该切片是否处于接触范围 if theta_local theta_in || theta_local theta_out continue; end k_slice1 calc_slice_stiffness(theta_local, rb(1), ra(1), E, nu, mn, z1, beta); k_slice2 calc_slice_stiffness(theta_local, rb(2), ra(2), E, nu, mn, z2, beta); inv_slice_pair 1/k_slice1 1/k_slice2; inv_k_pair inv_k_pair inv_slice_pair; end k_pair 1/inv_k_pair; % 多齿对并联叠加重复对相邻齿对执行上述过程并求和 Kmesh(j) k_pair; % 示例仅含单对齿 end这段代码里有两个关键判断。第一个是局部啮合角度的范围判断决定哪些切片真正承载。第二个是柔度合成方向先串联主动轮加从动轮再并联所有切片。有人会把切片方向先并联再串联结果会完全不同物理上也不对。实际的斜齿轮副在同一时刻会有多对齿同时参与啮合程序里需要在外层增加一个齿对循环把当前齿对和相邻齿对的贡献全部算完再累加。重合度由端面重合度和纵向重合度共同决定标准斜齿轮纵向重合度大于 1意味着至少有两对齿同时接触这个叠加不能省。4.4 计算效率的优化点如果切片数是 80扫描角度步长是 0.02 弧度一个完整啮合周期约 0.3 弧度那只有十几帧计算很快。如果扫描范围覆盖整个 2π帧数就变多单齿积分又是逐点做整体耗时可能达到几十秒到几分钟。一个有效优化是把齿廓积分点缓存起来。齿廓离散点只依赖齿轮几何不随啮合转角变化完全可以提前算好存成结构体数组。每次转角变化只更新积分上下限不必重复计算齿廓坐标。实测这个优化能把同一组参数的计算时间压到原来的三分之一以下。另外一个优化是切片刚度结果插值。预先在局部啮合转角维度上计算一组离散的切片刚度表再用一维插值获取任意转角的刚度。这样做可以在不损失精度的前提下大幅减少积分调用次数特别适合批量扫掠螺旋角和齿宽场景。5. 典型结果验证与常见翻车点程序写出来不等于算出来的东西能用。刚度模型最怕的是算完不知道自己算得对不对只知道画出来一条曲线。这一节给出一套能落地的验证思路和几个我踩过坑的典型问题。5.1 怎么验证程序算出来的值第一层验证是极限退化对比。把螺旋角设成 0程序应该退化成直齿齿轮模型。此时相位滞后项全部为零切片结果各片相同合成的刚度和等宽度直齿势能法结果应该一致。这一步能排查出切片法装配逻辑里的大部分错误尤其是相位偏移符号写反的问题。第二层验证是趋势对比。斜齿轮的时变啮合刚度曲线应该比直齿轮平坦波动幅度明显低于直齿轮。如果你算出的曲线在啮入和啮出位置出现大幅度锯齿多半是相位滞后没加或者切片合成方向搞错。第三层验证是数值量级。对于模数 2.5、齿宽 40mm 的钢制斜齿轮副单对齿刚度通常在每毫米几十万牛顿量级。算出来如果只有每毫米几万牛多半是积分单位或者几何量写错了。比如弹性模量用了 Pa 而长度用了 mm数值差千倍。这类错误在轮齿计算里非常常见我在代码里统一用 mm 和 MPa 为单位把弹性模量写成 206000让数值保持在好读的范围。5.2 切片法的网格敏感性切片厚度对结果的影响前面提过我再做一次更完整的量化描述。固定所有齿轮参数不变切片数量分别取 20、40、80、160观察总刚度在一个啮合周期内的最大变化量。20 片时曲线不光明显锯齿最大值和最小值的差比 160 片高出约百分之七。40 片已经把主要形状抓住但局部斜率仍有微小抖动。80 片以上基本完全重合。这是因为斜齿轮接触线是连续的斜线切片太粗相当于把连续接触线离散成阶梯状阶梯高度就是切片厚度乘上螺旋角的投影量。切片越细阶梯越接近真实斜线。这个误差不是随机误差所以不能用多次计算取平均来消除只能是加密网格。加密到一定程度后收敛此时结果可信。5.3 我在实际调试中遇到的四个典型问题第一个是主从动轮半径用混。基圆半径是各轮自己的相位滞后计算时也必须用对应轮的基圆半径而不是平均半径。斜齿轮副的主从动轮齿数不同基圆半径不同同一轴向位置上的角度滞后也不同混用后曲线左右不对称。第二个是压力角在法向和端面的转换。基圆半径计算必须用端面压力角而赫兹接触估算里的等效曲率半径用的是齿轮副在接触点的法向曲率。不少初学者把法向压力角当作端面压力角直接用结果基圆半径偏大刚度偏小。程序里我在派生参数区统一做好角度转换不把三角函数散落在各个函数里。第三个是齿根圆角半径给得不对。齿体弹性刚度对齿根圆角几何十分敏感圆角半径稍一变整体刚度就有几个百分点的变化。我建议直接按实际刀具参数或齿轮手册生成齿根过渡曲线不要用默认圆弧代替否则模型标定后用于其他齿数时会失真。第四个是切片厚度与转角步长不匹配。转角步长太粗会把啮入点的峰值削掉切片太粗会把斜齿轮的平滑特性抹掉。两者需要协调转角步长对应的接触点沿齿宽推进距离最好小于切片厚度否则相位滞后在离散化后被扭曲后处理时还会出现接触线断裂的假象。5.4 结果后处理的一个实用技巧算完一条刚度曲线先别看数值先画出来观察波形周期。正常斜齿轮副的刚度波动周期应该是基节对应的啮合周期也就是齿轮每转过一个齿距对应的时间。如果波形周期明显偏长或者出现拍频说明同时在计算多齿对时齿序处理有误。我习惯于把单对齿刚度和多对齿叠加后的总刚度画在同一张图里可以直观看出重合度的影响区域。单对齿刚度曲线的局部峰值位置和多对齿叠加曲线的叠加区应当对齐如果错开基本可以肯定是齿对索引或相位偏移的计算问题。6. 还能怎么扩展从单齿到全齿轮副和故障模型这套程序做到现在完成的是最基本的斜齿轮时变啮合刚度求解。真正的工程应用还往往会涉及几个扩展方向我把自己试过的路径整理一下给需要的读者一个参考。第一个扩展是加入齿廓修形。斜齿轮实际应用中几乎都有齿向修形或齿廓修形用来补偿受载变形和误差。修形量体现在齿廓离散点上也就是把理想渐开线沿法向偏移一定量。实现上只需要在生成齿廓点时做一个偏移计算不需要改动刚度积分框架。这个扩展对轮齿载荷分布分析特别有用可以快速评估不同修形量下的刚度波动幅值变化。第二个扩展是考虑齿轮误差和磨损。啮合刚度本身是几何量的函数齿距误差、齿厚偏差、齿面磨损都会改变单点的接触状态。在切片模型里可以给每一片单独设置齿距误差或磨损深度这样主程序扫描到该片时自动跳过或改变接触位置模拟实际运行状态下的啮合刚度变化。我用这个方式模拟过轻度磨损和齿面点蚀对刚度的削弱效果比整体降刚度系数丰富得多。第三个扩展是输出给动力学模型使用。时变啮合刚度通常作为集中参数模型里的一个周期性激励项程序最后可以把不同啮合转角的刚度值连同对应的时间轴一起导出成表格方便导入其他仿真工具。我在实际项目中会额外输出一个包含啮合位置、单对刚度、总刚度、重合度信息的日志文件后续做多体动力学分析时直接读取省掉二次手算。第四个扩展是把齿轮副的啮合相位做成并行计算。Matlab的并行计算工具箱可以直接把这个程序里的转角循环改成 parfor切片数量较大时加速效果明显。实测一个需要三小时完成的参数扫描任务改成并行后大约四十分钟跑完。不过要注意切片刚度表缓存逻辑在多线程环境下要处理好避免每个 worker 重复计算同一份齿廓数据。如果只是做算法验证算到第五节的位置就可以收工了。但要在实际工程里用我更建议把整个程序组织成一个函数输入齿轮参数和材料参数输出完整的刚度曲线数组。这样不管换齿轮模数还是换螺旋角都是一次调用的事后面接参数优化、故障仿真都有统一的接口。我自己跑工程案例时也是这么操作的前期多花一点封装时间后面省下大量重复劳动。