基于势能法与切片法的斜齿轮时变啮合刚度MATLAB计算

在过去的项目中，我多次被斜齿轮的时变啮合刚度（TVMS）问题卡住。它不像直齿轮那样，刚度的周期性变化简单清晰，斜齿轮因为螺旋角的存在，接触线是倾斜的，啮合过程是渐进发生的，刚度的变化曲线也显得更“圆润”。但如果真要去计算它，很多人第一反应是拉个三维有限元模型去“硬算”，可做参数研究和动态仿真优化时，每改一次齿宽或螺旋角都要重新建模，太不划算了。用MATLAB基于势能法加切片法来求解，既能保证精度，又能把计算成本压到一个普通工程工作站都能接受的范围，这也是我更推荐的方式。

我将整个计算过程拆成了三个层次：第一步，用势能法把单个轮齿的变形能折算成等效刚度；第二步，把斜齿轮沿齿宽切成若干薄片，用“局部看成直齿轮”的思想处理螺旋角带来的相位差；第三步，把所有切片和多对啮合齿的刚度按串并联关系汇总，得到以啮合相位为自变量的时变刚度曲线，再用傅里叶级数拟合出一个可直接用于动力学方程的表达式。

1. 不是玄学：“时变”这两个字到底来自哪里

很多人对“时变啮合刚度”的理解停留在“刚度随啮合位置变化”这个模糊概念上，但真正做动力学仿真时，这条曲线的形状直接决定了振动响应的频谱成分。啮合刚度的本质是轮齿抵抗弹性变形的能力，而齿轮在转动过程中，参与啮合的齿数和接触点的位置始终在变，所以综合刚度必然随时间周期变化。

先看直齿轮。直齿轮端面重合度一般在1.2到1.8之间，也就是说齿轮副的啮合过程会出现单齿对啮合区和双齿对啮合区交替。双齿啮合时，两条齿对分担载荷，综合刚度明显升高；单齿啮合时，载荷集中在一个齿对上，刚度和都变小。这种“高低不平”的阶梯变化正是齿轮啮合激励的主要来源，也是轮齿修形设计中要重点缓解的成分。

斜齿轮则不同。由于螺旋角使齿面在进入啮合时是斜着接触的，同一时刻在接触面上会有多个齿对参与，纵向重合度往往大于1。所以斜齿轮刚度曲线的特点是“起伏幅度比直齿轮平缓得多”，但仍然不是常数。更关键的是，接触线长度随啮合过程变化，这个长度直接影响赫兹接触刚度和单齿承载量，计算时必须把它还原出来。

再往细看，时变刚度的物理来源可以归纳为四块：轮齿弯曲变形、剪切变形、齿根局部弹性变形（有时归入基体柔度）、接触区域的局部挤压变形。其中接触变形对应赫兹刚度，弯曲和剪切变形对应齿轮本体的弹性储能，它们的变化规律各不相同，需要分开计算再组合。

所以在写求解程序之前，我的建议是先把目标明确：你是拿它做振动响应定性分析，还是想精确预测啮合激励力的幅值？如果只要定性规律，用简化公式就能满足；如果需要精确时变刚度作为非线性动力学模型的输入，那势能法加切片法是性价比最高的路线。我后面的程序也主要是按后者来设计的。

2. 势能法是怎么把轮齿拆成一根梁的

势能法的核心思想，是假设轮齿可以近似为一个变截面悬臂梁，齿根圆处固定，齿面载荷作用在渐开线齿廓的位置上。载荷在齿面上产生单位法向力，这个力使轮齿产生弯曲、剪切和轴向压缩变形，每一种变形都对应一项弹性应变能。

把应变能写出来之后，利用能量等效关系：等效刚度 k 等于 2倍的变形能除以载荷平方，也就是 W = F^2 / (2k)。把积分形式的应变能展开，就得到各刚度的倒数表达式。

2.1 弯曲刚度

弯曲变形是最主要的成分。假设载荷作用点到齿根的距离为 d，作用点处的齿厚为 h，齿廓上任意截面到载荷作用点的水平距离为 x，该截面的厚度为 h_x，截面惯性矩 I_x = (b·h_x^3)/12，其中 b 为齿宽。那么弯曲柔度可以写为：

1/k_b = ∫_0^d ( (x·cosα - h·sinα)^2 ) / ( E·I_x ) dx

式中的 α 是载荷作用点处的压力角，E 是弹性模量。这个积分要沿着齿廓方向逐步求解，所以程序里必须先把齿廓离散成一系列截面，每个截面都要算出到载荷点的距离和截面厚度。

这里有一个容易踩的坑：如果把齿根直接当成固定端，算出来的弯曲刚度会偏大，因为实际齿根是弹性圆角过渡，不是刚性夹持边界。工程上通常做法是保留过渡圆角区域的计算，同时把基体柔度单独补充进去。否则整条刚度曲线的数值会整体偏高，和有限元结果对不上。

2.2 剪切刚度和轴向压缩刚度

剪切变形主要发生在截面高度相对较大的区域，因为齿根处的剪应力集中比较明显。剪切柔度积分式是：

1/k_s = ∫_0^d ( 1.2·cos^2α ) / ( G·A_x ) dx

式中 G 是剪切弹性模量，A_x 是截面积。系数1.2是矩形截面的剪切修正系数，它是从抛物线剪应力分布推导出来的。如果做成其他截面形状，这个系数会有差异，但轮齿按矩形截面近似时，1.2是比较通用的取值。

轴向压缩柔度则对应齿面法向力在齿根方向的投影分量：

1/k_a = ∫_0^d ( sin^2α ) / ( E·A_x ) dx

很多人会忽略这一项，觉得压缩变形很小。实际对比下来，它的量级确实只有弯曲和剪切的几十分之一，但对刚度曲线形状还是略有影响。既然编写程序时积分路径是共用的，多算这一项成本很低，没有理由省略。

2.3 赫兹接触刚度

赫兹刚度代表齿面接触区的局部弹性趋近量，它不是柔度积分能解决的问题，需要按接触理论单独估算。对于线接触齿轮，使用著名的闭合解形式：

k_h = π·E·L / ( 4·(1 - ν^2) )

这里的 L 是接触线瞬时长度。这个式子假设接触区附近材料均匀、接触载荷按椭圆分布。实际接触载荷会影响赫兹接触宽度和趋近量，所以在更精细的模型里，k_h 不完全是个常数，而会随单位载荷轻微变化。不过大多数 TVMS 研究中，都把 k_h 近似为只依赖接触线长度和材料参数的量，这个简化在工程精度范围内是可接受的。

2.4 串联叠加规则

单个齿对沿法线方向的综合刚度，不能直接把这五部分相加，因为它们是沿变形传递路径“串联”的。整个齿对的变形等于弯曲、剪切、轴向压缩、赫兹接触和基体变形之和，所以柔度应该相加：

1/k_pair = 1/k_b + 1/k_s + 1/k_a + 1/k_h + 1/k_f

k_f 为基体柔度对应的刚度，反映了齿根过渡圆角以下弹性体对齿面位移的贡献。在不做有限元标定时，可以先用圆角区域延伸的简化近似公式估算，或者暂时忽略，但要明确知道这会带来多大的偏差。我自己的对比显示，忽略 k_f 可能使总刚度被高估6%到12%，具体大小与齿数、齿根圆角半径有关。

2.5 多齿对并联

设同一时刻有 n 对齿轮啮合，则总刚度是各齿对刚度的并联叠加：

k_total = Σ k_pair_i

因此在程序中，最外层的循环应该是“遍历当前啮合相位下所有参与啮合的齿对”，内层才是“沿齿廓的柔度积分”，或者切成薄片后的切片叠加。

3. 切片法：把斜齿“切直”再拼回去

直接求解斜齿轮的三维齿面弯曲问题，公式会变得非常复杂，因为轮齿不再是一个平面应变模型，螺旋角使得载荷位置沿轴向连续移动。切片法的思路很朴素：把齿宽 B 沿轴向切成 n 个薄片，每个薄片厚度为 Δz = B/n。由于薄片非常薄，每个切片上的几何关系可以近似看成一模一样的直齿轮剖面，只是它的“啮合时刻”相对整体参考面有一个轴向相位差。

3.1 切片相位差的计算

斜齿轮螺旋面的特点是：沿轴向移动单位距离，齿轮的端面齿廓相对于法向参考面转动一个角度。以基圆半径 r_b 为基准，轴向位置 z 处的相位偏移量可以写成：

δ(z) = z·tanβ / r_b

其中 β 是节圆上的螺旋角。这个偏移量的物理含义是，在同一个时刻，不同切片上的齿面并不在同一个端面投影位置，因此它们的渐开线展开状态是不同的。只要把这个偏移量叠加到齿轮的旋转角度上，就能把三维螺旋齿面当成一系列“不同初始相位”的直齿切片来处理。

我在程序里的实现方法是：先确定参考端面（z=0）上的啮合相位角 θ_0，然后对任意切片 z，实际用于几何计算的有效角为 θ_0 + δ(z)。注意这里的符号方向取决于螺旋角旋向，左旋和右旋要区分，否则切片刚度的分布会成镜像错位，导致接触线长度算错。

3.2 接触线长度的变化

斜齿轮的另一个关键几何量是接触线长度。直齿轮瞬时接触线平行于轴线，长度基本等于齿宽；斜齿轮的接触线是一条斜线，它从齿顶一端开始，朝齿根方向延伸，在啮合过程中长度先增大后减小。用切片法计算时，我实际上不需要专门去求接触线的解析表达式，只要在每个切片上判断“当前相位下，该切片是否处于啮合区域”，处于啮合区域的切片就贡献一个刚度，最后把所有贡献的切片刚度相加。这样，接触线长度被隐式地包含进来了。

例如，纵重合度 εβ 越大，同时处于啮合区的切片数就越多。当载荷从单侧进入时，参与切片数由0增加到 n，刚度平滑上升；离开啮合时再由 n 降到0。这个过程让斜齿轮的时变刚度曲线明显比直齿轮“圆润”。

3.3 切片和齿对的组合逻辑

我建议在程序中把啮合齿轮看成“主循环按齿对索引，内层按切片索引”。某一时刻可能存在两个甚至三个齿对同时啮合，那么每个齿对都要独立判断哪些切片处于啮合状态。齿对 i 的等效刚度为：

k_pair_i = Σ_j k_slice_ij

其中 j 遍历所有处于啮合状态的切片。得到所有齿对的等效刚度后，再用并联公式：

k_total = Σ_i k_pair_i

这种方式的优点在于代码结构和物理模型完全对应。之后如果想加入齿向修形，只需在每个切片刚度上乘以一个考虑修形量的权重系数，就能模拟不同切片的实际接触载荷和弹性变形分配，可扩展性非常好。

4. MATLAB 实现细节与核心代码段

搭建这个模型的MATLAB程序，并不需要特别高级的优化技术，但代码结构一定要清晰。整个程序可以分成四个模块：输入参数解析、齿廓离散与切片网格生成、单切片刚度积分、相位循环计算与拟合输出。

4.1 输入参数：先准备好这些量

注意，斜齿轮计算中很多几何量必须在端面坐标系里进行转换，比如端面压力角要写成：

alpha_t = atan( tan(alpha_n) / cos(beta) )

基圆半径是 r_b = 0.5·m_n·z / cos(alpha_t) 吗？这里要小心，端面模数 m_t = m_n / cos(beta)，所以基圆半径其实是 r_b = 0.5·m_t·z·cos(alpha_t) = 0.5·m_n·z·cos(alpha_t)/cos(beta)。这个转换错了，后面所有相位关系都会错。

4.2 齿廓离散和切片网格

我写程序时喜欢把每个切片的齿廓点存成结构体数组。对每个切片，渐开线的起点在基圆上，终点在齿顶圆附近，另外还要包含齿根过渡圆角。齿廓离散不光是几何上好看，它直接决定了柔度积分中的截面厚度 h_x 和截面积 A_x 能不能算准。

沿齿高的离散点数量建议不少于50个。对于斜齿轮，由于每个切片的有效相位不同，理论上每个切片都要独立离散一遍。如果齿数多、切片数也多，离散计算量会上升，所以更高效的办法是：只计算一个参考切片的离散几何参数表，其他切片通过角度差值复用同一套几何映射。这种方法只在变位系数和齿顶高系数不随轴向变化时成立，而正常齿轮设计就是这样，所以我实际代码里也采用了这个优化。

4.3 单切片刚度积分函数

下面的 MATLAB 核心函数用来计算某个切片上、某个啮合角位置的单齿对综合刚度。

function [k_slice, detail] = slice_stiffness(zeta, geo, mats) % zeta: 当前切片的有效啮合角 % geo: 包含齿数、模数、螺旋角、齿宽等结构体 % mats: 材料参数结构体 E = mats.E; G = mats.G; nu = mats.nu; b_slice = geo.B / geo.n_slice; % 单切片齿宽 alpha_t = atan(tan(geo.alpha_n)/cos(geo.beta)); r_b = 0.5 * geo.m_n * geo.z / cos(geo.beta) * cos(alpha_t); r_x = sqrt(r_b^2 + (r_b * zeta)^2); % 当前啮合点半径 alpha_x = atan2(r_b * zeta, r_b); % 该点压力角 % 预计算齿廓离散表 profiles：包含 x_c, h_c, I_c, A_c prof = geo.profiles; % n_dis x 4 矩阵 dmax = r_x - geo.r_root; % 载荷点到齿根距离近似 inv = zeros(geo.n_dis, 1); for i = 1:geo.n_dis x_c = prof(i,1); h_c = prof(i,2); I_c = b_slice * h_c^3 / 12; A_c = b_slice * h_c; inv_b = inv_b + ((x_c * cos(alpha_x) - h_c * sin(alpha_x))^2) ... / (E * I_c) * geo.dx; inv_s = inv_s + (1.2 * cos(alpha_x)^2) / (G * A_c) * geo.dx; inv_a = inv_a + (sin(alpha_x)^2) / (E * A_c) * geo.dx; end % 赫兹项，接触线长度取该切片有效宽度 kh = pi * E * b_slice / (4 * (1 - nu^2)); % 基体柔度由前述公式或外部函数计算 kf = foundation_stiffness(geo, mats, r_x); k_slice = 1 / (inv_b + inv_s + inv_a + 1/kh + 1/kf); detail = struct('kb',1/inv_b,'ks',1/inv_s,'ka',1/inv_a, ... 'kh',kh,'kf',kf); end

这个函数返回的就是“切成薄片后等效为一个直齿轮齿对”的刚度。真正用于整体刚度时，还要在外层叠加所有同啮合状态的切片，并进一步叠加多齿对并联。

4.4 主循环：相位扫描

外层主循环可以设置为从0扫描到2π，步长越细，曲线越光滑。但对于一个完整啮合周期，刚度以“一个齿距”为周期，因此更高效的写法是只扫描一个基节对应的角度区间，比如θ_span = 2π / z1。对斜齿轮，因为纵向重合度可能大于1，扫描区间要适当扩到两个到三个基节区间，以便观察多个齿对的交叠过程。我的经验是扫描角度取1.5个到2个基节，足够覆盖εγ不超过3的情况。

theta_cycle = linspace(0, 2*pi/geo.z1, N_step); K_total = zeros(N_step, 1); for it = 1:N_step theta0 = theta_cycle(it); K_sum = 0; for pair = 1:n_pair_max for j = 1:geo.n_slice zz = (j - 0.5) * geo.B / geo.n_slice; delta_phi = zz * tan(geo.beta) / r_b; % 相位偏移 zeta = theta0 + delta_phi; % 有效相位 % 判断是否处于啮合区间，可设一个区域上下限 if in_mesh_zone(zeta, geo) k_s = slice_stiffness(zeta, geo, mats); K_sum = K_sum + k_s; end end end K_total(it) = K_sum; end

这里要注意 pair 循环的边界设置。斜齿轮纵向重合度经常超过2，不能简单地和直齿轮一样只循环两对齿。合理的做法是根据总重合度εγ = εα + εβ确定最大可能同时啮合齿对数，再额外多留一对缓冲，避免漏算。

4.5 离散细节的收敛性

程序写完之后，第一件要做的事情是收敛性测试：把切片数从40加到160，观察刚度曲线的最大值和最小值变化。如果变化超过1%，说明切片数不够；如果变化很小，就选择那个“不再明显变化”的切片数。通常齿宽在30 mm，螺旋角15°左右时，80到120个切片已经足够。齿廓离散点也是类似的逻辑，50到100个点之间基本能收敛。

5. 把刚度曲线压缩成一个傅里叶公式

计算完成之后，你手里有一堆离散的刚度值：K_total(1)到K_total(N_step)。直接把这些离散数组丢给动力学软件也可以，但用起来很麻烦，特别是需要在Simulink中做连续仿真时，表查找反而容易引入插值误差。更好的做法是拟合出一个显式函数，把时变刚度变成光滑、可导的周期函数。

以啮合相位角 θ 为变量，刚度函数可以写成傅里叶级数：

K(θ) = K0 + Σ_{n=1}^{N} [ a_n·cos(n·θ) + b_n·sin(n·θ) ]

因为你本来就是在MATLAB里算的，直接调用 fft 或者 lsqcurvefit 都可以。我习惯先做 FFT，再截断到前8阶谐波。过低的阶数会保留太多波纹信息，过高的阶数会在数据边缘处带来“Runge”式振荡，8阶是一个比较稳妥的经验值。

下面的示例表展示了一次模拟计算的拟合系数，单位为N/μm。实际项目的系数务必用自己的几何参数重新计算，不能直接搬用。

拟合完成后，我会算三个指标：决定系数 R²、均方根误差 RMSE、最大局部误差。一般来说，R²都在0.995以上，最大局部误差控制在5%以内。对于齿轮动力学方程来说，这个精度完全足够，而且在刚度的一阶导数上都连续，不会在积分求解器中制造数值毛刺。

这个拟合公式的价值在于，后续做参数扫描时，不需要重新把整条曲线拿来查表，只需要更新十几个系数，就可以快速进入系统级仿真。如果你在开发齿轮传动系统的状态空间模型，这种表达方式也能直接参与矩阵组装，非常自然。

6. 验证与避坑实录

这个模型看着原理清楚，真正跑起来还是会遇到不少问题。我把个人项目里踩过的坑列出来，供大家对照排查。

6.1 刚度曲线出现“毛刺”或突变

如果曲线不是光滑周期变化，而是在某个位置突然掉下来或者跳上去，首先检查切片边缘的啮合判断是否过于生硬。切片在进入和退出啮合时，齿面接触是逐渐加载的，不是瞬间全接触。可以在切片边缘做线性过渡加权，比如引入一个0.05到0.15周期的渡越区，刚度乘以一个0到1的渐变系数，曲线会光滑很多。

6.2 赫兹刚度直接串联导致整体偏低

这是初学者最容易犯的错误。赫兹接触变形发生在接触表层，但它与齿轮体的弯曲变形并不是真正的“串联”关系。很多文献里为了工程简便还是按串联处理，这有一定合理性，但你需要知道这样算出来的总刚度会偏低。严格的做法是先按轮齿变形求出齿面法向位移，再把接触趋近量叠加上去，两者的变形耦合关系比简单串联更复杂。在工程精度要求不高时，串联近似法可以用，但对比实验数据时要留一个系数余量。

6.3 螺旋角旋向和相位方向搞反

切片相位项 δ = z·tanβ / r_b 的符号，取决于螺旋角旋向。左旋和右旋在同一个旋转方向下，齿面接触沿轴向的进给方向相反。如果程序里不区分，接触刚度的峰值位置会整个偏移一个切片相位，导致曲线形状完全错误。验证方法是：逐步增大螺旋角，观察刚度曲线从一个“直齿轮形状”平滑过渡到“圆润斜齿轮形状”，如果出现反向扭曲，多半是符号反了。

6.4 基体柔度对整体结果的量级影响

早期版本为了省时间，我没有单独算基体柔度，而是把齿根圆角直接当固定端。后来和实验结果对比，刚度偏高了8%左右，尤其在齿数少、齿根圆角小的齿轮上更明显。现在我的处理方式是，把齿根过渡圆角到外侧边界那一段位形单独建立二维弹性基础模型，用简单的平面应力公式估算基体柔度，再把结果作为常数项加进柔度求和里。这样对比验证结果吻合度好很多，程序运行时间也只增加了不到2%。

6.5 网格和点数不是越多越好

切片数从100加到200，结果几乎不变；但步长从0.01°缩到0.002°，计算时间涨了5倍，刚度变化只有0.3%。所以不要盲目追求网格密。判断收敛的方式很简单：“结果变化小于你使用的工程容差，就是够了”。对多数齿轮设计场景，总刚度变化1%以内的网格密度就可以接受。

6.6 如何验证程序正确性

我很建议用有限元做一个基准算例对比，不必做全齿三维网格，用一个单齿模型加周期性条件就能验证。比之前注意一点：有限元模型里轮齿的基体边界条件会影响刚度，两侧要是完全固定，刚度会比半轮体模型大不少。最好取一个大于两倍齿高厚度的齿圈壳体，再用对称约束处理，才能和切片法结果对齐。

另外，把螺旋角设成0°，你的斜齿轮程序应该退化为直齿轮程序，这可以作为一个强有力的“冒烟测试”。如果退化后结果和直齿轮经典公式对不上，问题大概率出在程序本身的切片逻辑上，而不是螺旋角处理上。

7. 这套模型还能往哪里延伸

算完基础的时变啮合刚度后，我发现这套框架的扩展能力比想象中强很多。第一个直接扩展方向是计及齿向修形和齿廓修形。切片法的天然优势就在于它已经按轴向切分了，只要在每个切片上引入一个修形量，用来修正该切片啮合位置的初始间隙，并通过间隙和载荷的迭代求出实际分担载荷，就能模拟修形后的刚度曲线变化。这个功能对于研究修形参数对传动误差和振动的影响非常有用。

第二个扩展方向是考虑变位系数的影响。改变变位系数会改变齿厚和齿根圆直径，也会改变啮合起始圆位置和压力角分布。切片法直接改输入几何参数即可，不需要重构任何算法逻辑。做齿轮传动方案优选时，可以批量扫参数，把整个设计空间的TVMS曲线都算出来，然后对比动态激励幅值。

第三个方向是和动力学方程耦合。拟合出的傅里叶系数可以直接组装到时变刚度微分方程里，后面再接一个模态求解器，就能分析系统的参数激励振动特性。比起从CAD模型导出再到有限元计算，这种纯数值路线在优化迭代场景下至少快一到两个数量级。

还有一点值得留意：如果齿轮转速很高，需要考虑齿面温升对弹性模量的影响，以及离心力引起的齿体变形。这些效应会把“静态”时变刚度变得与转速相关。切片法可以通过在刚度项中加入转速修正系数来近似处理，但这时模型参数必须依赖实验数据标定，不能纯靠公式推导。

一点实际体会

做这套模型的连续项目经验告诉我，时变啮合刚度的价值不在于单个数值算得多准，而在于它能否反映不同设计变量下刚度的相对变化趋势。用MATLAB把势能法和切片法结合以后，参数扫描变得很轻松，每次修改设计都能快速看到刚度曲线形态的改变。如果你也想沿着这个方向做，建议第一版程序先保持朴素，把单个切片的刚度计算验证准确，再加傅里叶拟合，最后再去做修形和变位的扩展。咬住“先直齿后斜齿、先单体后整对”的节奏，后面的路会顺很多。