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 重复计算同一份齿廓数据。

如果只是做算法验证，算到第五节的位置就可以收工了。但要在实际工程里用，我更建议把整个程序组织成一个函数，输入齿轮参数和材料参数，输出完整的刚度曲线数组。这样不管换齿轮模数还是换螺旋角，都是一次调用的事，后面接参数优化、故障仿真都有统一的接口。我自己跑工程案例时也是这么操作的，前期多花一点封装时间，后面省下大量重复劳动。