行星齿轮箱的动力学仿真，第一步永远是刚度。无论是做固有特性分析、动态响应计算还是齿面载荷分配，时变啮合刚度都是方程里躲不开的核心参数。这次要用程序解决的就是行星传动中最典型的一对啮合：行星轮与内齿圈构成的啮合副，也就是标题里的"行星齿轮内啮合齿轮副"。
我写的这套程序采用势能法，定位在健康齿状态，也就是齿面没有裂纹、没有剥落、没有过度磨损，只考虑标准渐开线齿形下的理论啮合刚度。程序比较大一部分精力花在了"精确渐开线齿形"上，不是拿圆弧拟合那种偷懒做法，而是老老实实把渐开线方程代进几何求解，从根到顶的齿廓坐标都按渐开线规律生成。这对后续算弯曲刚度、剪切刚度、齿基柔度都至关重要，齿形一旦是近似值，刚度曲线的形状和峰值就会跟着飘。
文章后面我会把程序的核心思路、内啮合特有的几何坑、五个刚度分量的计算流程、多齿啮合叠加逻辑全部摊开讲，也会把调试中遇到的几个典型问题列成速查表。适合正在做齿轮动力学、轮齿修形或者故障诊断相关课题的人参考，尤其是准备从外啮合刚度程序转向内啮合副计算、又不太想一上来就啃有限元的新手。
1. 程序定位与整体设计思路
1.1 这个程序到底解决什么问题
行星齿轮副的时变啮合刚度，本质上就是一对齿在从进入啮合到退出啮合的过程中，单位齿宽上的法向载荷与法向变形之间的比值。这个比值不是常数，因为啮合点沿着啮合线移动，轮齿悬臂梁的有效长度在变，接触点处的曲率半径也在变，重合度又决定了同一时刻有几对齿同时参与承载。所以刚度是啮合相位的函数，画出来是一条周期性曲线。
把这个刚度曲线做准，是后续一切动力学的前提。系统矩阵里的啮合刚度项直接决定固有频率和振型；时变刚度作为参数激励，又是齿轮副啮合振动响应和边频带形成的根源。很多人在做行星齿轮箱故障诊断时，拿健康状态下的时变啮合刚度基线去对照损伤状态，那个"健康齿"基线，来源就是这套理论计算程序。
内啮合副和外啮合副的区别在于，内齿圈的齿顶圆在齿根圆的内侧，轮齿朝中心方向伸出，几何上整个就反过来了。直接套外啮合的公式，齿顶圆半径、基圆半径、啮合线长度这些符号全乱。程序里所有几何量都按内啮合重新推导过，包括重合度和啮合点的坐标，这样才能保证行星轮-齿圈这对副的刚度算出来是正常量级。
1.2 为什么选势能法，而不是有限元
齿轮啮合刚度主流算法有三类：有限元法、解析公式法、势能法。有限元最准，能做齿根过渡圆角处的应力集中，能考虑齿圈轮缘的整体弹性，问题是要建精细网格、设接触对、逐相位求解，一个啮合周期算几十个位置，模型量大、耗时长，做参数扫描或者后续叠加载荷谱就不太现实。解析公式法快，但大多基于经验系数修正，齿形一变位、齿数一变，误差立刻放大。
势能法属于半解析方法。它把轮齿当成变截面悬臂梁，把齿基当成弹性体，把接触区当成赫兹接触，逐个计算各个变形分量对应的刚度，再用柔度串联的方式叠加。这个思路物理意义清楚，计算量又小，适合做工程研究。更重要的是，势能法天然跟"渐开线几何"深度绑定——悬臂梁积分要沿齿高方向逐段取截面，每个截面的厚度必须由渐开线方程算出来；载荷作用点的压力角也要由渐开线展角确定。所以势能法程序写得好不好，很大程度上取决于渐开线齿形建模得细不细。
我最终选择了势能法，目标就是把计算速度和几何保真度平衡起来。一个啮合周期离散成两三百个相位点，每个相位点算两到三对齿，MATLAB单循环跑完也就几十秒，换参数重算几乎零成本。这对课题组做参数化分析来说太重要了。
2. 精确渐开线齿形与内啮合几何建模
2.1 渐开线方程怎么代进程序
渐开线的本质是一条直线绕基圆纯滚动时，直线端点的轨迹。写成参数方程就是：
x = r_b * (cos(t) + t * sin(t)) y = r_b * (sin(t) - t * cos(t))这个 t 是滚动角，它的正切值等于接触点到基圆的切线段长度除以基圆半径。程序中所有关键量——接触点半径、载荷角、任意截面的齿厚——都要先从这条方程推出来。比如齿面上某一点对应的压力角 α_k，满足：
cos(α_k) = r_b / r_p其中 r_p 是该点的极径。知道极径之后，这个点在渐开线上的展开角就是 inv α_k = tan(α_k) - α_k。
精确齿形的含义就是：给定齿数、模数、压力角、变位系数，先算出基圆半径和分度圆弧齿厚，然后用渐开线方程逐点生成齿面坐标，从齿根工作圆一直扫到齿顶。这样齿面上任何一点的半径、齿厚弧长、压力角都是真实值，不是拿圆弧替代。
齿厚计算同样不能省。外齿轮某一半径处的弧齿厚是：
s_r = r * [s_ref / r_ref + 2 * (inv α_ref - inv α_r)]s_ref 是分度圆上的弧齿厚，r_ref 是分度圆半径，α_r 是半径为 r 处的压力角。内齿轮的表达式形式相同，但齿厚在径向方向上从齿根向齿顶收缩，符号处理必须按内齿轮的"齿顶在内、齿根在外"重新排列。程序里专门写了一个齿厚求解子函数，输入半径输出弧齿厚，悬臂梁积分时直接调用。
2.2 内啮合副的几何反直觉点
内啮合副的几何，新手第一次做容易栽跟头的地方是齿顶圆和齿根圆的半径关系。以模数2、齿数60的内齿圈为例，按标准直齿轮：
- 齿顶圆半径 ra = m * (Z - 2) / 2 = 58 mm
- 齿根圆半径 rf = m * (Z + 2.5) / 2 = 62.5 mm
- 基圆半径 rb = m * Z * cos(20°) / 2 ≈ 56.38 mm
看看这三者的关系：齿根圆在外侧最大，齿顶圆在内侧最小，基圆介于两者之间。这跟外齿轮"齿顶圆最大、齿根圆最小"的印象完全相反。所以内啮合的啮合线长度公式必须单独推导，不能照搬外啮合。
我程序里实际啮合线长度用的是内啮合专用表达式：
g_alpha = sqrt(ra_in^2 - rb_in^2) - sqrt(ra_p^2 - rb_p^2) + a * sin(α)其中 ra_in、rb_in 是内齿圈参数，ra_p、rb_p 是行星轮参数，a 是中心距。这个式子的符号很关键：内齿圈那一项带正号，行星轮那一项带负号，中心距投影项也带正号，三个量叠加得到实际啮合线长度。如果符号搞反，算出来的重合度要么小于1，要么出现负数，一看就知道几何没理顺。
2.3 重合度与啮合区间的确定
重合度公式是实际啮合线长度除以基圆齿距：
ε = g_alpha / (π * m * cos α)拿模数2、行星轮20齿、内齿圈60齿、压力角20度的组合来说，实际啮合线长度算出来大约是15.86 mm，基圆齿距5.90 mm，重合度约2.69。这意味着传动过程中大部分时间有两对齿同时啮合，中间一小段区域会达到三对齿同时承担载荷。
这个重合度对刚度曲线的影响非常直接。重合度大于2之后，刚度曲线不再是经典的"双齿区-单齿区-双齿区"的V形或U形，而变成"两齿区-三齿区-两齿区"的交替平台。内齿圈齿数越多，重合度越大，曲线波动率反而越小。很多文献里说内啮合行星级刚度波动小于外啮合，根源就在这里，不是材料变好了，是几何的重叠效应把刚度拉平了。
程序里做啮合区间划分时，就是按照基圆齿距为基本步长，把一个啮合周期（基节长度）离散成若干相位点。对每个相位点，先判断当前有几对齿处于啮合线内：把第 i 对齿的进入点位置和退出点位置求出来，再看接触点坐标落在哪些齿对的啮合区间里。这一步用逻辑判断实现，循环体很轻，效率很高。
3. 势能法计算时变啮合刚度的核心实现
3.1 五个刚度分量的物理来源和公式
势能法把单对齿的啮合柔度拆成五个部分：赫兹接触柔度、弯曲柔度、剪切柔度、轴向压缩柔度、齿基柔度。每一部分的物理来源不同，最后通过柔度相加再取倒数得到单齿对啮合刚度。
赫兹接触刚度描述的是两齿面在接触点处的局部弹性压陷。对钢制齿轮，使用经典的线接触赫兹公式：
1/k_h = 4 * (1 - ν^2) / (π * E * L)这里 E 是弹性模量，ν 是泊松比，L 是齿宽。公式里不含载荷大小，因为线接触刚度的线性化结果就是这个形式，实际齿面载荷分布不均带来的修正可以留到后续做修形分析时再细化。
弯曲刚度是核心，它把轮齿当成一个固定端在齿根、自由端在接触点的变截面悬臂梁。接触点处的法向载荷 F 分解成垂直于齿中心线的分量 F * cos(α_k) 和平行于中心线的分量 F * sin(α_k)。沿齿高方向从齿根积分到接触点：
1/k_b = ∫[ (d - x) * cos(α_k) - h_x * sin(α_k) ]^2 / (E * I_x) dxd 是接触点到齿根的径向距离，x 是积分位置的径向坐标，h_x 是载荷作用线到该截面的偏心距，I_x 是该截面的惯性矩。I_x 由齿厚决定：I_x = L * s_x^3 / 12，其中 s_x 是半径 x 处的弧齿厚。注意这里截面惯性矩按矩形截面近似，齿宽 L 为厚度方向，齿厚方向为弯曲方向。这个近似在齿数较多时很准，齿数太少、齿形矮胖时需要再斟酌。
剪切刚度用截面平均剪应力近似：
1/k_s = ∫[1.2 * cos^2(α_k)] / (G * A_x) dx系数1.2是矩形截面剪应力分布不均匀的修正因子，G 是剪切模量，A_x = L * s_x 是截面积。轴向压缩刚度对应平行于齿中心线的载荷分量：
1/k_a = ∫[sin^2(α_k)] / (E * A_x) dx这一项通常比弯曲刚度小一两个数量级，但既然做理论程序，就一起算上，免得凑不齐总柔度。
齿基柔度描述的是齿根以下轮体弹性变形引起的那部分附加柔度。我用的是Sainsot经验公式体系，把齿根过渡区几何无量纲化，再查经验系数。表达式形式比较长：
1/k_f = cos^2(α_k) / (E * L) * [L_*(u_f/S_f)^2 + M_*(u_f/S_f) + P_*(1 + Q_*tan^2(α_k))]其中的 L_、M_、P_、Q_是多项式插值系数，u_f 是载荷作用线到齿根圆角危险截面的距离，S_f 是该截面处的齿厚。内齿轮的 u_f 计算和外齿轮不一样，因为内齿轮的危险截面在半径更大的齿根侧，载荷作用点在内侧，两者的相对位置符号是反的。程序里这里单独写了条件分支，专门处理内齿圈的齿基柔度计算。
3.2 单齿对柔度组合与多齿叠加逻辑
单对齿的总柔度，等于两个轮齿各自柔度之和再加接触柔度：
1/k_pair = (1/k_bp + 1/k_sp + 1/k_ap + 1/k_fp) + (1/k_br + 1/k_sr + 1/k_ar + 1/k_fr) + 1/k_h然后取倒数得到单齿对啮合刚度。这里的下标 p 代表行星轮，r 代表内齿圈。为什么是柔度相加而不是刚度相加？因为变形在载荷方向上串联叠加，力相同、变形相加，柔度满足加法关系，直观类比成两根弹簧首尾相接。
同一时刻若有 N 对齿同时啮合，它们分担的是同一个法向载荷，但每对齿的变形各自独立，刚度上并联。所以系统总刚度为所有参与啮合齿对的刚度之和：
k_m = Σ k_pair_i这里隐含一个假设：各齿对均匀分担载荷。工程上这个近似在健康齿、未修形条件下是可接受的，考虑齿面误差和载荷分配不均时需要引入传递误差协调方程，程序暂时没有走到这一步。
多齿叠加的实现并不复杂，关键是相位对齐。啮合周期内，第 i 对齿比第 i+1 对齿晚进入啮合一个基节距离。程序以行星轮的角位置为自变量，把每对齿的接触点沿啮合线的位置写成：
s_i(t) = s_entry_i + (t - t_entry_i) * v_relativev_relative 是啮合点沿啮合线的移动速度，等于基圆线速度。对每个离散相位 t，只需要判断每一对齿的 s_i 是否落在 [0, g_alpha] 区间内，落在区间内的齿对计入并联求和。
3.3 主程序流程与关键代码逻辑
完整程序分七步，这里把流程主线列出来：
1 输入基本参数：模数m、行星轮齿数Z_p、内齿圈齿数Z_r、压力角α、变位系数x_p、x_r、齿宽L、材料参数E、ν 2 计算几何参数：分度圆、基圆、齿顶圆、齿根圆、中心距 3 计算实际啮合线长度g_alpha、基圆齿距pb、重合度ε 4 生成渐开线齿廓坐标：行星轮齿面、内齿圈齿面 5 将啮合线[0, g_alpha]离散为N个相位点 6 对每个相位点调用单齿对刚度子函数，得到N组单齿对刚度 7 按重合度区间判断参与啮合的齿对序号，并联叠加得到时变啮合刚度曲线核心的单齿对刚度子函数，内部核心逻辑是：
% 输入：接触点半径 r_k、压力角 alpha_k d = r_k - r_root; % 齿根到接触点的径向距离 x = linspace(0, d, 200); % 沿齿高积分网格 sum_b = 0; sum_s = 0; sum_a = 0; for i = 1:length(x) r_x = r_root + x(i); % 当前截面半径 s_x = toothThickness(r_x); % 渐开线齿厚函数 I_x = L * s_x^3 / 12; A_x = L * s_x; h_x = 载荷线偏心距(r_x, alpha_k); sum_b = sum_b + ((d - x(i))*cos(alpha_k) - h_x*sin(alpha_k))^2 / (E*I_x) * (x(2)-x(1)); sum_s = sum_s + 1.2*cos(alpha_k)^2 / (G*A_x) * (x(2)-x(1)); sum_a = sum_a + sin(alpha_k)^2 / (E*A_x) * (x(2)-x(1)); end k_b_inv = sum_b; % 齿基柔度、赫兹柔度分别计算后相加代码里特意把积分网格取200个点，这个数量经过实测足够保证收敛。网格太少，刚度曲线在啮合区边缘会出现锯齿；网格再多，计算时间线性增长但精度改善可以忽略。齿厚函数 toothThickness 是渐开线几何模块的核心，内部用 inv 函数换算压力角与齿厚的对应关系，内齿圈调用时传入一个方向参数，保证齿厚沿径向收缩方向正确。
4. 调试记录与常见问题排查
4.1 刚度值整体偏大或偏小
程序跑出来的第一版曲线，如果平均刚度比文献值高出一截，优先检查齿基柔度有没有漏算。很多人第一次写势能法程序，弯曲、剪切、轴向三项都齐了，唯独把齿基柔度当成小量忽略，结果总刚度偏大10%到15%，峰峰值形状看着没问题，绝对值对不上，对不上就说明少了一路柔度并联。
反过来，如果刚度整体偏小，检查齿根圆半径取值。外啮合程序里齿根圆半径用的是分度圆减齿根高，内齿圈的齿根圆半径却是分度圆加齿根高加顶隙。齿根圆取错，悬臂梁长度 d 就变长，弯曲柔度积分区间拉大，刚度就掉了。这个符号问题非常隐蔽，我调试时靠打印中间变量才发现，一行一行对比 d 和 r_root 的数值，才定位到是齿根半径的正负号取反了。
另外一个常见原因是积分网格里的齿厚函数没有按渐开线算。如果齿厚直接取分度圆弧齿厚常数，等于把变截面梁当等截面梁算，弯曲刚度会明显偏差。齿数越多齿厚变化越平缓，误差越小；齿数少、齿形相对较矮时，这个误差会把刚度曲线抬高20%以上。
| 现象 | 可能原因 | 排查方法 |
|---|---|---|
| 刚度整体偏大 | 齿基柔度缺失或系数取错 | 单独输出各柔度分量，核对 k_f 占比 |
| 刚度整体偏小 | 内齿圈齿根圆半径符号取反 | 打印 r_root，与标准公式对照 |
| 曲线台阶明显 | 积分网格过粗 | 网格加到300点，观察曲线是否收敛 |
| 峰值点抖动量异常 | 载荷作用点压力角用了名义值 | 确认 α_k 由接触点半径实时求解 |
4.2 曲线不连续或跳变
时变啮合刚度曲线最常见的问题是齿对交替处出现台阶跳变，比如某几个相位点刚度突降，形成肉眼可见的凹口。这个凹口往往不是物理现象，而是啮合齿数判定逻辑出错。
我用的判定方法是：把每对齿的接触点位置 s_i(t) 算出来，判断是否在 [0, g_alpha] 范围内。问题出在边界处理上：接触点正好等于0或等于 g_alpha 的那一个相位点，前后两个程序分支一个算进去了，一个没算进去，导致齿对数突然少了一对或突然多了一对。解决方式很简单，把区间判定改成闭区间包含端点，并且在端点处做一个加权处理：接触点位于啮合线端点时，该齿对刚进入或刚退出，载荷为零，刚度贡献也趋向零，所以无论算不算入总刚度结果都应连续。我最终把边界处接触点到端点距离小于一个微量的齿对强制按零刚度处理，曲线立刻变顺了。
还有一个坑是接触点沿啮合线的移动速度。如果程序里把行星轮角度增量换算成弧长时，用的半径是分度圆半径而不是基圆半径，位移量和实际啮合线长度就对不上。验证方法很简单：相邻两个相位点的接触点间距，理论上应当等于该相位区间内的基圆弧长，而不是分度圆弧长。这个错误会在刚度曲线的周期性上暴露出来——相邻两个齿距的曲线形状看起来相似但又不完全一致，波长对不齐。
4.3 与有限元结果对标
程序写完后，建议拿一组简单参数跟有限元结果对一下。我自己用ABAQUS算过模数2、齿数比20比60、齿宽20mm的单齿对刚度，对比结果显示势能法程序误差大概在5%左右，曲线趋势完全一致。误差主要来源是齿基柔度，有限元算出的齿根附近变形比经验公式偏大，因为经验公式的外推系数最早是依据外齿轮标定的。
对标时注意一个细节：有限元模型要固定内齿圈外圆节点，模拟行星齿轮箱中齿圈与箱体连接的工况。如果不约束齿圈外圆，或者把内齿圈当成自由环处理，轮缘整体变形会把齿基柔度放大，刚度偏低。这个约束条件的选择对对标结果影响很大，比网格密度的影响还明显。程序里对应的是齿基柔度公式中 u_f 的取值，也就是从齿根圆到齿圈外圆的径向距离，外圆半径越大，齿圈轮缘越厚，齿基柔度越小。
注意：有限元对标时，内齿圈外圆约束方式必须和程序中的齿基柔度外推假设一致，否则误差大得离谱，很可能不是程序写错而是边界条件没对上。
5. 实操心得与扩展方向
5.1 健康齿基线的用处和边界
健康齿状态下的时变啮合刚度曲线算出来后，最直接的用处就是作为故障诊断的基线。比如齿根裂纹，裂纹一来会导致齿根截面的有效惯性矩下降，弯曲刚度局部减小，反映到刚度曲线上就是在对应裂纹位置的啮合相位出现局部凹坑。把健康齿基线存成基准数组，再做带损伤齿的势能法程序，两个结果相减，就能得到刚度降的时域分布，定位到是哪根齿、哪个啮合区段出了问题。
不过要提醒一点，健康齿基线对几何参数非常敏感。程序里如果输入了变位系数，齿厚的实际分布已经偏离标准齿轮，刚度曲线会明显改变。所以建立基线时参数一定要拿实际齿轮图纸来填，不能随便套默认值。加工误差、修缘、齿向鼓形这些实际状态，在健康齿程序里统统没考虑，它们本质上已经属于"非健康"状态了。
5.2 程序后续可以怎么扩展
这套程序的扩展空间还比较大。一条路是加均载系数：把多齿并联的简单求和改成变形协调方程，考虑各齿对的传递误差协调分配载荷，这样更接近真实工况。另一条路是加齿根过渡曲线：现在我用的是理论齿根圆角，实际滚齿或插齿的过渡曲线形状不同，对齿根应力集中和齿基柔度都有影响，把刀具齿顶圆角代入过渡曲线方程，精度还能再往上提。再就是加裂纹模型，目前业界常用的做法是把裂纹等效为截面惯性矩折减，在悬臂梁积分中按裂纹深度和位置修正 I_x，这部分我后面也在做，等程序稳定了再单独写一篇详细过程。
5.3 内齿轮建模的几条经验总结
最后说几条个人经验，都是踩过坑之后记下来的。
第一，内啮合程序的所有几何公式，先画一张从内齿圈中心出发的简图再写代码。齿顶圆在内、齿根圆在外这个反直觉关系，看公式一百遍不如把图画一遍。程序里凡是涉及半径比较的地方，统统加上注释说明大小关系，防止未来某天回来改代码时又被符号绕晕。
第二，渐开线函数的弧度制问题。inv(α) = tan(α) - α，这个式子里的 α 必须用弧度。程序里输入压力角时用度数，后面所有计算转弧度，但到了 inv 函数内部，很容易把 trans 函数和输入单位混在一起。我一开始就栽在这里，一度刚度曲线歪得完全不像样。现在所有角度量在函数入口统一转弧度，命名上直接带 rad 后缀，彻底杜绝这类问题。
第三，齿基柔度公式里的经验系数，首次使用不要直接信任默认参数。我的做法是先拿一对简单的外啮合副（比如20齿对20齿）跑一个基准算例，跟文献里的结果对比，确认整个势能法框架没问题，再切换到内啮合工况。这样如果内啮合结果异常，问题基本可以锁定在内啮合专有的几何处理部分，排查范围小很多。这套思路让我的调试效率提升了不少，如果你也要从零写一套齿轮刚度程序，建议照这个顺序走。
