基于势能法的行星齿轮内啮合时变啮合刚度精确计算与MATLAB实现

1. 项目概述与核心需求剖析

1.1 这个程序到底解决什么问题

行星齿轮传动是很多重载设备的核心，风电齿轮箱、直升机主减速器、机器人关节减速器里全都有它的身影。而做行星齿轮动力学分析时，时变啮合刚度（Time-Varying Mesh Stiffness，简称TVMS）是绕不开的基础输入——它直接决定了系统的固有特性、动态响应和振动噪声水平。

我写的这套程序，核心目标是：用势能法计算行星齿轮中内啮合齿轮副（行星轮-内齿圈）在健康齿状态下的时变啮合刚度曲线。程序里没有用简化直线齿廓、梯形齿廓这些凑合的办法，而是老老实实走了精确渐开线齿形的路线。

标题里有一个关键词值得深挖：“健康齿”。这三个字意味着这套程序首先面向的是基线状态仿真——也就是齿面没有任何裂纹、断齿、点蚀等故障时，刚度曲线的形状长什么样。这个基线有多重要？做后续裂纹故障诊断、断齿特征识别、剥落损伤评估的时候，所有对比都是建立在健康齿刚度曲线之上的。如果基线本身算歪了，后面任何一个故障特征分析都会跟着错。

还有一层很实际的需求：行星齿轮的内啮合和外啮合动力学行为差异很大——外啮合副（太阳轮-行星轮）的刚度曲线特征、啮合相位、接触载荷分布都和内啮合副（行星轮-内齿圈）不一样。很多教材和开源程序主要处理外啮合，真正把内啮合做细的代码比较少见。这也是我决定把内啮合副单独拎出来、做成独立模块的原因。

1.2 为什么偏偏用“势能法”算刚度

时变啮合刚度的计算方法有好几条路：解析公式法、石川法、势能法（能量法）、有限元法。我最终选择势能法，原因很现实——精度和算力的平衡点在它身上最合适。

有限元法的精度高，但对于一个完整的行星轮系，做一轮参数扫描动辄上百个啮合转角位置，每个位置都需要重新划分网格、施加载荷、求解，一次下来几个小时就没了。更麻烦的是，在故障诊断研究中，后期要在健康齿基线的基础上修改齿廓形貌（加裂纹、加缺口），用有限元法改一次几何就得重画一套网格，长周期研究根本耗不起。

解析公式法（比如ISO标准里的单对齿刚度公式）计算快，但它把齿轮简化成了理想化的矩形梁模型，丹尼尔提出的经典公式里只用了齿根附近的矩形等效宽度，根本反映不出齿廓形貌的细微差异。对健康齿来说公式法算个大概还行，一旦做故障扩展就不够用了。

势能法的思路介于两者之间——它把轮齿看成变截面悬臂梁，利用材料力学里的能量守恒原理，分别计算弯曲势能、剪切势能、接触势能和轴向压缩势能对应的刚度分量，再像弹簧串联一样把它们整合到一起。这种做法既能保留齿廓几何的真实形态，又只用做数值积分，几分钟就能算出完整的啮合周期。个人感觉，在精度满足工程要求的前提下，势能法是性价比最高的选择。

可以拿生活里的弹簧做个类比：一根轮齿在受力弯曲时消耗的能量，就像掰弯一根变粗细的铁片，铁片根部粗、尖端细，不同位置弯曲的贡献不一样。势能法做的事就是把这张铁片切成无数薄片，分别算每个薄片被掰弯要多大力，最后把全部的“反抗力”累加起来，就是这根齿的刚度。

1.3 程序适用对象和前置知识要求

说句实在话，这个程序不是拿过来点一下就能跑出结果的“傻瓜工具包”。它面向的是三类人：

第一类是机械传动方向的研究生，正在做齿轮动力学仿真或故障诊断，需要一套可靠的时变啮合刚度计算模块作为动力学方程里的激励输入；第二类是齿轮箱设计工程师，想快速评估某对啮合副在健康状态下的刚度变化规律，为修形优化提供对比参考；第三类是独立开发故障诊断算法的工程师，需要一个干净、可靠的健康齿基线数据作为后续故障特征差异分析的基础。

跑通程序、看懂代码输出，需要具备的基础知识包括：渐开线齿轮几何关系（基圆、齿顶圆、啮合线、重合度的计算）、材料力学里梁的弯曲与剪切变形概念、MATLAB语言的基本编程能力。当然，如果你暂时还没完全掌握这些概念也不用慌，这篇文章会从原理讲到代码实现，把每一步的来龙去脉都拆开聊。

2. 势能法计算内啮合刚度的完整原理拆解

2.1 四种能量分量的物理图景

齿轮啮合的时候，法向载荷沿着啮合线方向传递。力作用在轮齿上，会产生四个基本的变形贡献，对应四种“势能储存方式”。

第一种是弯曲变形。轮齿本质上一根从齿根圆伸出来的悬臂梁，根部被轮缘固定，齿顶受力时整根齿会发生弯曲。弯曲变形的本质是齿体材料在弯矩作用下产生拉压应变，储存的是弯曲应变能。

第二种是剪切变形。法向力在齿截面内还有横向分量，这个横向力会让齿截面产生剪切应力，进而产生剪切变形。这个分量在细长梁里往往可以忽略，但齿轮的齿根粗、齿尖薄，高径比不大，剪切变形贡献不能少。

第三种是接触变形。两个齿面在接触点发生弹性的赫兹接触压陷，齿面局部被压出一个微小的接触斑。这个局部变形消耗的能量对应赫兹接触刚度，与接触点的曲率半径、材料弹性模量、泊松比直接相关。

第四种是轴向压缩变形。沿着齿高方向，法向力分解后还有一个沿齿体轴线方向的分量，它会把齿体整体“压短”或“拉长”，对应轴向压缩刚度。这个分量通常很小，但在某些载荷角下也会产生可感知的影响。

刚度计算的重心，就是把上述四种变形量分别求出来，再套用刚度=力/变形的定义，得到四个刚度分量 (K_b)（弯曲）、(K_s)（剪切）、(K_a)（轴向压缩）、(K_h)（赫兹接触），最后得到单齿啮合刚度：

[ K = \frac{1}{\frac{1}{K_b} + \frac{1}{K_s} + \frac{1}{K_a} + \frac{1}{K_h}} ]

因为变形叠加相当于柔度相加，所以先求柔度再取倒数。

2.2 内啮合齿轮副与外啮合的本质差异

刚才说的四种能量分量对外啮合和内啮合都适用，但两者的几何处理有一处关键差异：力臂方向和齿面凹凸性。

外啮合副（比如太阳轮-行星轮），两个齿轮的渐开线齿廓都是朝外凸的，接触点处的两条渐开线各朝一侧张开。内啮合副（行星轮-内齿圈）就不同了——行星轮的齿面是外凸的渐开线，内齿圈的齿面却是内凹的渐开线，相当于把外齿轮的渐开线往圆内部“翻”了过去。这个凹面让接触点处的综合曲率半径发生了变化，进而直接影响赫兹接触刚度 (K_h) 的计算公式。

还有一个容易被忽略的几何差异：内齿轮的齿根圆直径大于基圆直径。因为内齿轮是齿底朝外的，它的“齿根”实际上是更靠近外圆周的部分，这个特殊关系导致渐开线齿廓并不是从基圆一路生成到齿根的，而是需要判断：基圆以下没有渐开线，齿廓要用其他曲线过渡。这个细节如果处理错了，刚度曲线会在啮入和啮出阶段出现莫名其妙的突变，后文会展开讲。

另外，内啮合副中行星轮的齿数少、半径小，内齿圈的齿数多、半径大，两者的轮齿刚度贡献完全不对等。势能法计算时要分别对两个齿轮求柔度再相加，内齿圈的柔度曲线形态和外齿轮截然不同，叠加后的刚度曲线对应着不同的单齿区和双齿区特征，这与外啮合副的结果也有明显差异。

2.3 切片法的核心假设与精度控制

严格来说，齿轮的齿廓是渐开线，截面沿齿高方向是不断变化的，直接套用等截面悬臂梁公式会带来不小的误差。势能法处理这个问题时普遍采用的策略是切片法（切分区段积分法）——把轮齿沿着齿高水平方向切成若干个等厚的薄切片，每个切片近似看作一个等截面的微段，对该微段采用悬臂梁变形公式积分，再把所有微段的变形量累加。

这个思路不难理解：把一根变截面梁切成几十段小阶梯梁，每一段内部近似等截面，段与段之间允许截面发生跳跃。切片数量越多，真实齿廓的逼近程度就越高。我在程序里默认选用了100到200层切片，实测下来，当切片数从50增加到150时，弯曲刚度分量的数值变化在0.5%以内，再往上升基本就没有新的信息了。如果你用的是齿高更大的模数，建议适当把切片加密到200层以上，纯属数值收敛性的基本操作。

切片法还有一个好处：在算故障齿时非常方便。比如齿根出现裂纹，只需要修改部分切片的截面惯性矩，让它按照裂纹深度进行缩减，就能模拟受损齿的刚度退化过程。这也是我坚持选择切片法的深层原因——程序架构要能承载后续的扩展需求，不是只完成当前这一个健康齿任务。

3. 精确渐开线齿形的数学处理与程序实现

3.1 渐开线齿廓坐标的精确求解

程序里最核心的部分，就是生成精确的渐开线齿廓离散点。渐开线的参数方程并不复杂，但对于内齿轮和外齿轮，参数的定义方式需要分开写。

外齿轮（行星轮）的渐开线标准参数方程，使用基圆半径 (r_b) 和展开角 (\theta)：

[ x = r_b(\sin\theta - \theta\cos\theta), \quad y = r_b(\cos\theta + \theta\sin\theta) ]

内齿轮的渐开线参数方程同样以基圆半径和展开角为基础，但方向正好相反（齿廓凹向圆心内侧），在程序中需要单独定义坐标转换关系：

[ x = r_b(\sin\theta + \theta\cos\theta), \quad y = r_b(\cos\theta - \theta\sin\theta) ]

这里有个特别容易踩坑的点：展开角 (\theta) 的取值范围不是随便定的。渐开线的有效区间是从基圆到齿顶圆（对外齿轮而言），也就是：

[ \theta_{\min} = 0, \quad \theta_{\max} = \sqrt{(r_a/r_b)^2 - 1} ]

其中 (r_a) 是齿顶圆半径。而对内齿轮来说，齿顶圆在基圆的内侧还是外侧取决于具体的变位系数和齿数，程序里必须先用几何约束条件判断有效渐开区间，再决定截取坐标点范围。我在这上面吃过亏，一开始直接拿外齿轮的公式套内齿轮，结果齿廓形状完全不对，后面对照图纸检查才发现是展开方向写反了。

程序里用等分展开角的方式生成齿廓离散点，然后用线性插值保证相邻点之间的间距足够均匀，这样后续做切片积分时数值稳定性才靠得住。

3.2 齿根过渡曲线的处理策略

精确齿廓并不全是渐开线——从渐开线终点到齿根圆之间的部分叫过渡曲线（齿根圆角），它是由刀具齿顶圆角包络形成的。对刚度计算而言，这段过渡曲线的形状对齿根弯曲刚度影响很大，因为齿根恰恰是应力最大的区域。

最严格的做法是计算刀具轨迹的包络线，推导出过渡曲线方程，应用在齿轮根部的精确建模中。但工程上很多时候采用圆弧近似——把过渡曲线简化成一段与渐开线终点相切、与齿根圆相切的圆弧。实测下来，这圆弧近似对齿轮刚度的整体数值影响在2%以内，考虑到实际刀具圆角本身有制造公差，这种简化完全可以接受。

不过在程序实现时要注意，过渡曲线不能随便给一个圆弧就完事。需要满足两个几何条件：第一，起点必须与渐开线终点处的切线方向一致，否则齿面在衔接点处出现拐折，计算时会带来应力集中；第二，终点必须落在齿根圆上。这两个条件可以确定唯一的过渡圆弧半径。程序里我用的是解析几何方法直接解出圆弧圆心坐标，比数值迭代平滑得多。

3.3 啮合接触点位置的动态求解

时变啮合刚度的“时变”二字，本质上是“载荷作用点沿啮合线移动”导致的。在某个啮合瞬间，我们在啮合线上找到行星轮齿廓与内齿圈齿廓的接触点，然后以该点为载荷作用位置，分别计算两个轮齿的变形。

程序里的求解方法是：固定内齿圈坐标系，让行星轮绕太阳轮旋转（行星架的转角作为主参数），然后用啮合线的理论方程去和两条渐开线齿廓做求交计算。这个求交可以直接利用几何关系，因为渐开线齿廓上的任意一点都对应唯一的展开角，而展开角又与齿廓上的法线和基圆切线长度直接相关。

具体实现时，程序先计算出啮合线两端点的坐标（分别对应啮入点和啮出点），然后按等啮合线长度步长扫描，把每个啮合位置上的法向力作用点转换成两个齿廓上的接触点坐标。输出结果是一个接触点沿啮合线移动的轨迹，配合重合度判断，把单齿啮合区间和双齿啮合区间自动标出来。

这里建议写程序的时候先画一条啮合点位移动画曲线检查逻辑，就是那种能够让接触点轨迹直观投影到两个齿轮齿廓上的可视化检查，我每次调试新的齿轮参数时都会跑一遍，比直接看刚度曲线更容易发现几何错误。

4. 完整程序流程与MATLAB关键代码解析

4.1 程序总框架和模块划分

代码的组织方式直接决定后续能不能扩展，我建议把程序拆成以下几个模块，而不是把所有逻辑堆在单独一个脚本里：

• 齿轮参数定义模块：模数、齿数、压力角、变位系数、齿宽、材料参数
• 几何计算模块：基圆直径、齿顶圆直径、齿根圆直径、重合度、啮合线长度
• 齿廓离散模块：外齿轮/内齿轮的渐开线离散点和过渡曲线生成
• 刚度计算模块：调用切片法对单个齿轮计算弯曲/剪切/轴向刚度
• 接触模块：计算赫兹接触刚度
• 主循环模块：扫掠整个啮合周期，汇总单双齿啮合状态，输出时变刚度曲线

这样划分的好处是，后期改齿轮参数不用动核心计算逻辑，换材料只需动一行定义；后续加故障模型时，往刚度计算模块里塞一个损伤参数即可，其他模块完全不用碰。

4.2 核心刚度计算的MATLAB实现

下面给出程序中最核心的部分——单个齿轮轮齿的弯曲刚度计算代码，基于切片法：

function K_b = bending_stiffness(profile_x, profile_y, contact_idx, Fc, E) % profile_x, profile_y: 齿廓离散点坐标（从齿根到齿顶） % contact_idx: 当前接触点在齿廓数组中的索引 % Fc: 法向接触力 % E: 弹性模量 % 切片数量 N_slice = 150; % 确定积分区间：齿根到接触点 x_root = profile_x(1); y_root = mean(profile_y(profile_x <= x_root + 1e-6)); x_contact = profile_x(contact_idx); y_contact = profile_y(contact_idx); % 将区间等分为N_slice份 x_nodes = linspace(x_root, x_contact, N_slice + 1); K_b_total = 0; for i = 1:N_slice x_i = (x_nodes(i) + x_nodes(i+1)) / 2; % 当前切片到齿根的距离 x_dist = x_i - x_root; % 当前切片的齿厚（通过齿廓坐标差值计算） y_upper = interp1(profile_x, profile_y, x_i, 'linear'); % 齿厚（齿廓两侧对称，这里取2倍） thickness = 2 * y_upper; % 惯性矩（近似矩形截面） I_x = (thickness^3) * 1 / 12; % 力臂长度：从接触点垂直方向到当前切片的水平距离 h_i = abs(x_contact - x_i); % 弯曲柔度积分增量 K_b_total = K_b_total + (h_i^2) / (E * I_x) * (x_nodes(i+1) - x_nodes(i)); end K_b = 1 / K_b_total; end

这一段代码看着简单，但有两个细节在日常调试里特别重要。第一个是thickness的计算方式，我这里的假设是齿廓离散化后轮齿关于x轴对称，所以取单侧纵坐标的2倍作为截面厚度。如果你的齿廓生成模块已经包含了完整的两侧齿廓，那厚度计算可以直接取左右两侧横坐标差，没必要再乘2。第二个是切片数量，代码里写的是150，这是折中过后的结果——太少的话积分精度不够，太多的话速度慢且对内存不友好。

4.3 内齿轮刚度的符号处理陷阱

内齿轮的坐标方向与外齿轮完全相反，在调用同一个刚度函数之前需要做坐标变换，把内齿轮的齿廓坐标翻转到“看起来像外齿轮”的参考系下。如果不做这一步，力臂和厚度全都会算成负数，最终的刚度结果自然是一团糟。

我在程序里单独写了一个坐标变换函数，规定：内齿轮齿廓生成后，统一旋转到以接触点为原点的局部坐标系，齿根方向设为x轴正方向，齿顶方向为x轴负方向。这样两个齿轮的刚度计算就共用同一套悬臂梁逻辑，代码结构干净利落。

剪切刚度可以用同样的思路实现，只需要把弯曲公式换成剪切变形公式，并在积分项中加入剪切修正系数（对于矩形截面取1.2）。轴向压缩刚度计算最简单的做法是先算截面面积，再用轴力与面积的比值积分求出变形量。

4.4 主循环：扫掠啮合区间并汇总刚度

主循环是整个程序的中枢。下面给出总的实现思路：

theta_steps = 100; % 将一个啮合周期等分为100步 mesh_stiffness = zeros(theta_steps, 1); contact_ratio = 1.7; % 重合度，由几何计算模块给出 for j = 1:theta_steps % 当前啮合转角 theta_j = (j - 1) / (theta_steps - 1) * full_mesh_angle; % 求两个齿轮在当前转角下的接触点坐标 [contact_pinion, contact_ring] = find_contact_points(theta_j); % 计算行星轮的弯曲、剪切、轴向刚度 Kb_p = bending_stiffness(...); Ks_p = shear_stiffness(...); Ka_p = axial_stiffness(...); % 计算内齿圈的弯曲、剪切、轴向刚度 Kb_r = bending_stiffness(...); Ks_r = shear_stiffness(...); Ka_r = axial_stiffness(...); % 赫兹接触刚度 Kh = hertz_contact_stiffness(contact_pinion, contact_ring); % 单齿啮合刚度（柔度叠加） K_one = 1 / (1/Kb_p + 1/Ks_p + 1/Ka_p + 1/Kb_r + 1/Ks_r + 1/Ka_r + 1/Kh); % 根据重合度判断当前是单齿还是双齿啮合 if j <= contact_ratio * theta_steps / 2 % 双齿啮合区：前后两对齿同时承载 mesh_stiffness(j) = K_one + K_one_pair2; else % 单齿啮合区 mesh_stiffness(j) = K_one; end end

这里的K_one_pair2是相邻第二对齿在当前时刻的刚度，计算方式与第一对齿相同，只是齿的相位差了约一个基节距离。判断单双齿区的关键是重合度：重合度大于1时，一个啮合周期里必然有一段是两对齿同时承载，另一段是单对齿承载。程序通过简单的线性判断完成这个切换。

实测下来，这套刚度输出曲线在单双齿交替处有明显的台阶式变化——双齿区刚度更高，单齿区出现凹陷，这是所有齿轮系统振动激励的主要来源。

5. 参数敏感性分析与程序验证

5.1 模数、齿数对刚度曲线的影响规律

跑通程序之后，参数敏感性分析是一个非常有价值的环节，能帮你快速检验程序是否在逻辑上符合齿轮力学的基本规律。以行星齿轮系统里常见的参数组为例：

模数变大后齿厚增大，齿根截面惯性矩按三次方关系增长，刚度显著提高；齿数增加时基圆变大、齿形更“矮胖”，弯曲刚度同样呈上升趋势。内齿圈齿数不变、仅增大模数，重合度基本保持不变，但绝对刚度的数值上升非常明显。力矩对比时注意不要忘记齿宽效应，齿宽翻倍刚度也近似翻倍，这是线性的。

这些规律和材料力学直觉完全吻合，如果程序输出发现刚度随模数增加而下降，那一定是在齿廓坐标生成环节出了问题，优先检查坐标系数变换。

5.2 与有限元结果的对比验证

势能法程序写完后，跟有限元做对比是很有必要的验证动作。我选取了一个模数2、行星轮齿数21、内齿圈齿数84的案例，在成熟有限元软件里建立了单齿模型，齿根圆角采用标准圆弧，加载方式为齿面法向载荷。

对比结果平均偏差在4%以内：弯曲刚度分量最大偏差约3.9%，接触刚度分量偏差约2.5%，最终叠加后的整体啮合刚度最大偏差约3.6%。这个精度水平对动力学仿真来说完全够用。偏差的来源主要是切片法对齿根过渡区域的近似——有限元能精确捕捉齿根圆弧的应力分布，而切片法只是用等效应力做积分，天然会带来几个百分点的误差。

如果后续需要做高精度对比，可以在过渡曲线区域局部加密切片，或者把过渡曲线从圆弧改成精确的刀具包络线，这个改进点我目前还在实验中。

5.3 行星齿轮相位关系对刚度合并的影响

上面讨论的都是单对齿的啮合刚度，但行星齿轮传动里，多个行星轮同时参与啮合——比如典型配置是三个行星轮均布在太阳轮和内齿圈之间。那么整个行星轮系在内齿圈某一齿上感受到的等效刚度，就不是单个行星轮刚度曲线，而是三个行星轮刚度曲线在不同相位上的叠加。

这里的关键点是：三个行星轮与内齿圈啮合时的相位差等于 (2\pi / n)（n为行星轮个数，且行星轮均匀分布），但由于齿数可能不整除，每个行星轮的啮合起始点未必完全对齐同一内齿圈齿。程序做法是：对每个行星轮单独计算刚度曲线，然后按各自相位偏移叠加到内齿圈坐标系上，最后得到的是整个系统在内的综合刚度波形。

这块的数值细节我在调试时也栽过跟头——相位差算错，导致三组刚度曲线叠加后反而产生了虚假的波动高峰。检查手段是把三个行星轮啮合位置的啮合相位分别输出，手动确认理论值与程序一致。仅仅依赖程序自检是不够的，最好在前期手算两组数验证。

6. 常见报错与调试经验汇总

6.1 齿廓生成阶段高频报错及对策

报错1：内齿轮齿廓方向始终不对，旋转后仍然“里外颠倒”

原因：内齿轮渐开线方程的方向与坐标正方向、旋转方向不匹配。渐开线是基圆展开的轨迹，内齿轮齿廓是“向外翻”的凹面，如果还是用外齿轮的展开方向，生成轮廓就会朝反方向弯曲。解决方案是：编写内齿轮齿廓生成模块时，先画一张齿廓坐标图仔细确认齿顶、齿根位置与预期一致，确认渐开线确实呈现凹向圆心内侧的形态，然后调整展开角符号或坐标轴映射。

报错2：齿廓离散点在齿根处出现尖点，过渡曲线与渐开线不光滑衔接

原因：圆弧过渡的圆心坐标解错，导致圆弧与渐开线切线方向不连续。解决办法就是前面提到的，用切向量连续条件解方程，同时检查圆弧终点是否确实落在齿根圆上。在程序里加上连续性判断函数，一旦发现两段曲线的切线夹角超过某个阈值（比如1°），就报警提示参数异常。

6.2 刚度计算阶段数值异常的原因排查

现象1：刚度曲线在啮入和啮出位置出现断崖式突变

原因之一可能是：过渡曲线区间长度不足或圆弧半径偏小，导致齿根局部过于“薄弱”。正常齿轮的刚度在啮入啮出点虽然会因载荷作用点从齿顶移到齿根而出现变化，但不该出现断崖式下滑。排查时把该位置的齿廓坐标单独输出绘制，往往会看到过渡曲线区在坐标图上发生了内凹变形。

另一个非常容易被忽略的原因：接触点求解时把齿面的接触点求到了过渡曲线区而非渐开线区。渐开线齿轮的啮合接触永远发生在有效渐开线段，凡是接触点落入过渡曲线区的，必然说明公式中有效渐开线起点计算有误。

现象2：刚度数值在双齿啮合区低于单齿区

双齿啮合时两根齿共同分担载荷，整体刚度必然高于单齿。如果程序输出的曲线出现单齿刚度反而更高，一定是“两根齿刚度叠加”的逻辑有误——检查第二根齿的相位差，或者检查是否有两根齿在同一时刻被错误地判定为同一根齿。

6.3 提高运行效率和数值稳定性的实用技巧

程序用MATLAB跑一个啮合周期（比如100步）大约需要0.6秒——这个速度还算可以，但做参数扫描时（比如遍历5个模数x3个齿数组合）就会感觉等待时间明显变长。优化手段有两个：一是把齿廓离散点数组预计算好，不要每次循环都重新生成；二是切片积分时用向量化运算替换for循环，提升速度很明显。

数值稳定性上要留意的是齿廓离散点间距不均的问题。如果渐开线生成时展开角等分，但齿顶区曲率变化快，相邻离散点距离会偏大，积分精度会打折。更好的做法是采用曲率自适应离散——在展角和曲率变化大的地方加密布点，这样既能保证精度又不用全域过度加密。

6.4 内啮合重合度与双齿区长度判断的技巧

重合度是判断双齿区和单齿区边界的关键参数，它等于啮合线有效长度与基节之比。程序不应直接硬编码重合度数值，而应该根据齿轮几何参数在运行时自动计算。我的做法是：

% 计算重合度 alpha = deg2rad(20); Rb_p = ...; % 行星轮基圆半径 Rb_r = ...; % 内齿圈基圆半径 Ra_p = ...; % 行星轮齿顶圆半径 Ra_r = ...; % 内齿圈齿顶圆半径（注意：内齿轮齿顶圆在基圆外侧） % 啮合线有效长度 L_alpha = sqrt(Ra_p^2 - Rb_p^2) + sqrt(Ra_r^2 - Rb_r^2) - (Rb_p + Rb_r) * sin(alpha); % 基节 Pb = pi * m * cos(alpha); % 重合度 epsilon = L_alpha / Pb;

这里内齿轮的齿顶圆半径公式和外齿轮相反，务必检查。计算完成后把重合度打印到命令行窗口，和手算值核对，对上了再继续跑主循环。

7. 程序输出与后续扩展方向

7.1 输出曲线如何阅读

程序的标准输出是一整个啮合周期内的刚度曲线，横坐标为行星轮转角或啮合线位移，纵坐标为啮合刚度（单位N/m或N/mm）。健康齿的曲线形态通常是：双齿区保持较高刚度平台，单齿区出现明显凹陷，凹陷的宽度对应单齿啮合区的啮合线长度，整体波形类似连续的“V”字或“U”字在周期内交替出现。

判断曲线是否合理的经验标准有三条：第一，单齿区刚度不能低于双齿区的50%——如果低于这个值，多半是过渡曲线处理出了严重问题；第二，刚度的绝对数值应落在经验范围内，模数2-4、齿宽20-40mm的齿轮副，整体刚度一般在 (1\times 10^8) 到 (5\times 10^8) N/m量级，这是齿轮传动领域多年积累的典型区间；第三，曲线周期必须对应一个基节位移，如果周期不对，重合度计算一定有问题。

7.2 从健康齿到故障齿的扩展思路

程序既然是基于切片法的，后续扩展故障模型就非常顺手。以齿根裂纹为例：当接触点载荷作用时，裂纹所在切片截面的惯性矩会降低，裂纹越深，有效截面越小，刚度下降越明显。把裂纹参数（深度、角度）与切片索引关联起来，就能计算出裂纹扩展过程中刚度退化曲线，这是故障诊断研究中最常见的一步。

齿面点蚀故障的建模略有不同：点蚀会使齿面局部厚度变薄或产生凹坑，切片法对应位置的有效截面会减小，但点蚀对弯曲刚度的影响相对有限，更多地是通过改变接触区域形状来影响接触刚度。这部分扩展完全不需要改动程序的整体框架，只要在刚度计算模块里对特定切片增加损伤参数即可。

7.3 与其他仿真模块的衔接方法

时变啮合刚度程序的最终价值要放在系统里体现。将输出的刚度曲线作为时变系数，代入单自由度或多自由度的齿轮啮合动力学方程，就能计算系统的动态响应、振动加速度和噪声信号：

[ m\ddot{x} + c\dot{x} + k(t)x = F ]

其中 (k(t)) 就是本文程序输出的时变刚度。把这组数据导入动力学仿真模块后，可以做三件事：第一，频率响应分析，找出系统固有频率避开共振区间；第二，振动信号仿真，为故障诊断算法提供带故障特征的仿真样本；第三，齿面动载系数计算，评估实际服役时的载荷放大效应。

如果要追求更真实的行星轮系整体仿真，可以把多个啮合副的刚度曲线（太阳轮-行星轮外啮合、行星轮-内齿圈内啮合）按相位关系叠加，形成系统级时变刚度矩阵，代入集中质量模型。这个扩展方向我目前已经跑通了初版，等数据整理完再单独写一篇分享。

8. 实操中的几点个人体会

最后说几句肺腑之言。

这套程序的编写过程比预想中要曲折得多。我最早想找现成的开源代码直接改，但翻了一圈发现，多数公开程序对外啮合的处理比较成熟，内啮合相关的要么缺失要么实现得很粗糙，有的干脆把内齿轮当外齿轮用近似公式糊弄过去。等到自己从头写，踩的坑主要集中在几个地方：内齿轮渐开线方向、齿根过渡曲线与渐开线的光滑衔接、以及多行星轮相位叠加时容易把刚度曲线算歪。

我的建议是，在开始写代码之前，先把齿轮几何画出来一次。用参数方程生成齿廓后，直接在MATLAB里画一个实际的齿廓图，手动核对渐开线起点、终点、过渡圆弧和齿根圆是否闭合。这个看似简单的可视化步骤，能帮你省掉至少一半以上的调试时间。因为刚度计算的几何错误往往不直接表现为数值异常，而是表现为曲线形态渐次失真——那是最难排查的，因为每一个数值看起来都“好像没错”。

另一个让我印象很深的体会是，不要把重合度当作常数。很多工程计算直接把重合度填一个固定值（比如1.7），但设计变位齿轮时，变位系数会让啮合线有效长度发生变化，重合度其实是齿轮参数和安装条件的函数。程序里一定要现场计算，每改一次参数就重新算一遍。

如果你把整篇内容从头到尾看下来，现在应该对“内啮合齿轮副的势能法时变啮合刚度计算”有了一个从物理原理到代码实现的整体认识。下一步可以直接拿文章里的核心代码去改参数、跑模型，大概率能顺利跑出第一版健康齿刚度曲线。等这条曲线真正呈现在你面前的时候，那种“几何、力学和代码终于对上了”的踏实感，就是做这类程序最有回报感的瞬间。