行星齿轮箱的动力学仿真第一步永远是刚度。无论是做固有特性分析、动态响应计算还是齿面载荷分配时变啮合刚度都是方程里躲不开的核心参数。这次要用程序解决的就是行星传动中最典型的一对啮合行星轮与内齿圈构成的啮合副也就是标题里的行星齿轮内啮合齿轮副。我写的这套程序采用势能法定位在健康齿状态也就是齿面没有裂纹、没有剥落、没有过度磨损只考虑标准渐开线齿形下的理论啮合刚度。程序比较大一部分精力花在了精确渐开线齿形上不是拿圆弧拟合那种偷懒做法而是老老实实把渐开线方程代进几何求解从根到顶的齿廓坐标都按渐开线规律生成。这对后续算弯曲刚度、剪切刚度、齿基柔度都至关重要齿形一旦是近似值刚度曲线的形状和峰值就会跟着飘。文章后面我会把程序的核心思路、内啮合特有的几何坑、五个刚度分量的计算流程、多齿啮合叠加逻辑全部摊开讲也会把调试中遇到的几个典型问题列成速查表。适合正在做齿轮动力学、轮齿修形或者故障诊断相关课题的人参考尤其是准备从外啮合刚度程序转向内啮合副计算、又不太想一上来就啃有限元的新手。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 对齿比第 i1 对齿晚进入啮合一个基节距离。程序以行星轮的角位置为自变量把每对齿的接触点沿啮合线的位置写成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齿跑一个基准算例跟文献里的结果对比确认整个势能法框架没问题再切换到内啮合工况。这样如果内啮合结果异常问题基本可以锁定在内啮合专有的几何处理部分排查范围小很多。这套思路让我的调试效率提升了不少如果你也要从零写一套齿轮刚度程序建议照这个顺序走。