企业官网建设流程全解析

1. 这个题目到底在算啥：健康齿轮的时变啮合刚度

做齿轮箱振动分析这几年，最容易被新手忽略的是啮合刚度其实不是一个常数。齿轮一转，参与啮合的轮齿对数和接触点位置都在变，啮合副的等效刚度自然也在跟着变。基于传统材料力学势能法的健康齿轮时变啮合刚度数值分析，就是把“健康齿轮”这条基线先用解析加数值的手段算清楚：不考虑裂纹、点蚀、断齿这些故障，只按完整轮齿的几何形状和弹性变形，得到一条随啮合周期稳定波动的刚度曲线。

为什么这条曲线这么重要？齿轮动态特性方程里，时变啮合刚度是核心的周期系数。它直接决定齿根动载荷、轮体振动、噪声辐射，也决定系统会不会出现参数共振。故障诊断里更是把健康齿轮的刚度曲线当参考基准，实际齿轮一有裂纹，曲线会在特定啮合位置出现明显的局部凹坑，再往后的边带、谐波分析都建立在这条基线上。所以健康状态的刚度不是“算个大概就行”，而是要尽量算得准、算得可复现。

适合谁参考？刚接手齿轮动力学课题的研究生、做减速机或新能源电驱系统NVH的工程师、以及想自己搭故障诊断仿真平台的开发者，都可以从这套方法入手。它不像有限元那样需要完整三维模型和复杂网格，用材料力学和数值积分就能得到与有限元非常接近的结果，而且计算速度快到可以做成批量扫参工具。

2. 势能法：把轮齿拆成一根变截面悬臂梁

2.1 为什么选势能法而不是直接上有限元

有人会问：现在有限元这么成熟，为什么还要用传统材料力学？我的经验是，有限元做一次精确的啮合刚度确实不难，难的是“变参数扫描”。齿轮模数、齿数、变位系数、齿宽一换，模型就要重建，或者至少重新划分接触区域网格；要模拟一个完整啮合周期内不同接触位置，还要在多个相位下做静力接触分析，模型规模大、求解时间长、后处理也繁琐。

势能法的优势在于把轮齿抽象成一根固定端在齿根、自由端承受接触力的变截面悬臂梁，再分别计算弯曲、剪切、轴向压缩和赫兹接触四种变形所储存的应变能，通过能量与刚度之间的关系反推出等效啮合刚度。整个过程只依赖几个几何参数，换一组参数就是改几行代码的事。而且它把变形来源拆得清清楚楚，能直接看出弯曲刚度占多大比例，这对后续分析裂纹位置的影响非常有利。

这套方法的局限也要心里有数：它基于平面假设和梁理论，对齿根过渡曲线、轮体柔度的处理比较粗。齿数很少、齿宽很大或轮辐结构很特殊的齿轮，纯势能法会低估齿根附近变形，这时候需要额外引入轮体柔度修正项，或者用有限元标定系数。对常规直齿轮副而言，直接算已经够用。

2.2 四个能量分量与刚度公式

设齿面法向载荷为F，轮齿总的弹性变形由四部分叠加：赫兹接触变形、弯曲变形、剪切变形、轴向压缩变形。因为有变形能U和刚度k的关系U = F²/2k，柔度可以写成1/k = 2U/F²。每个变形机制彼此串联，总柔度等于各项柔度之和：

1/k_m = 1/k_h + 1/k_b1 + 1/k_s1 + 1/k_a1 + 1/k_b2 + 1/k_s2 + 1/k_a2

下标1、2分别代表主动轮轮齿和从动轮轮齿。如果考虑轮体弹性，还要加上1/k_f1和1/k_f2两项。

赫兹接触刚度反映接触点处局部弹性压陷，像两个曲面在载荷下的挤压变形。对直齿轮副，接触线沿整个齿宽方向，近似有：

k_h = π E b / [4(1 − ν²)]

这里E是弹性模量，b是有效齿宽，ν是泊松比。如果两个齿轮材料不同，公式里的E取等效弹性模量E*，也就是两个齿轮材料参数的调和形式：1/E* = (1−ν₁²)/E₁ + (1−ν₂²)/E₂，再代入时系数要相应调整。为了方便，大多数同材料钢制齿轮副直接取E = 2.06×10⁵ MPa。

弯曲刚度是重点。把齿根截面固定，接触力作用在距离齿根为d的位置，任意截面x处的弯矩为F cos α₁ (d − x)，其中α₁是载荷方向与齿对称中心线的夹角。于是：

1/k_b = ∫₀ᵈ cos²α₁ (d − x)² / (E I(x)) dx

I(x)是该截面绕中性轴的惯性矩。直齿轮某截面若厚度为2h(x)，则I(x) = 2b h(x)³/3。

剪切刚度对应剪力引起的切向变形：

1/k_s = ∫₀ᵈ 1.2 cos²α₁ / (G A(x)) dx

系数1.2是矩形截面剪切形状因子，G为剪切模量，A(x) = 2b h(x)为截面面积。

轴向压缩刚度对应载荷沿齿高方向的分量：

1/k_a = ∫₀ᵈ sin²α₁ / (E A(x)) dx

这三组积分是势能法的核心。α₁的严格取值应该随接触位置变化，因为接触点沿齿廓移动时，法向力相对轮齿中心线的角度会变。工程实践中很多人直接取分度圆压力角α = 20°替代，cos²20°约等于0.88，影响不算太大。但对结果要求高时还是应该做逐点几何插值，我后面会讲到怎么处理。

2.3 轮体柔度要不要算

悬臂梁模型只把齿根截面以上当变形体，实际上轮体也会有弹性变形。齿根受到弯矩后，轮缘和轮毂部分会发生微量扭转变形，这个量对刚度贡献在短粗齿上尤其明显。经典处理方法是在总柔度中追加一项轮体柔度1/k_f，常见的有Sainsot给出的多项式经验公式，还有ISO齿轮承载能力标准里的参数化方法。

我自己的判断标准是：当齿数大于30、齿宽与模数比在8到20之间时，轮体柔度占比通常控制在10%以内，去掉影响不大。但如果是少齿数、大模数、薄轮缘这类结构，轮体柔度可能占到总量的20%以上，这时候必须加修正项。对本文的健康齿轮算例，我选择先不展开轮体柔度，专注把轮齿部分积分算准；需要更严苛对比时再补上。

3. 几何建模：从齿轮参数到可积分的齿廓

3.1 关键参数与渐开线方程

标准外啮合直齿轮的几何参数很固定：模数m、齿数z、分度圆压力角α_p、齿顶高系数h_a*、顶隙系数c*。分度圆半径r_p = m z/2，基圆半径r_b = r_p cos α_p，齿顶圆半径r_a = r_p + h_a* m，齿根圆半径r_f = r_p − (h_a* + c*) m。这些半径决定了积分边界。

渐开线齿廓是核心。任意半径r_K处的压力角α_K满足cos α_K = r_b/r_K，对应的渐开线函数为inv α_K = tan α_K − α_K。在半径r_K处，轮齿的弧齿厚可以由分度圆齿厚按渐开线展角换算得到：

s(r_K) = r_K [s_p/r_p + inv α_p − inv α_K]

其中s_p = π m/2是标准齿轮分度圆上的弧齿厚。这个式子的物理含义是：半径变化后，渐开线轮廓让齿厚按展角差重新分布。有了s(r_K)，任意截面齿厚度就是s(r_K)，半厚度h(x) = s(r_K)/2。

为什么强调渐开线函数？因为整个势能法对齿厚分布很敏感，尤其是靠近齿根处的截面惯性矩。如果把齿形简单当成梯形，齿根刚度会明显偏大，算出来的TVMS曲线会和真实齿轮差出一截。

3.2 齿根圆和基圆的相对位置是个坑

很多初学者最容易在这里出错。齿数少的时候，基圆半径反而大于齿根圆半径，也就是r_b > r_f，这意味着从齿根圆到基圆这一段并不存在真正的渐开线，而是齿根过渡曲线，通常由刀具圆角切出来。此时齿廓的渐开线部分是从基圆才开始向上的。

渐进线积分时，如果直接套用渐开线齿厚公式算到r_f以下，会把不存在的几何硬算出来，结果虽然数值连续，但物理上不对。我处理这类情况有两种办法：第一种是只从r_b积分到载荷点，齿根到基圆的过渡段单独用直线近似；第二种是编程时做一个保护判断，当r < r_b时固定取齿根附近的最小齿厚，避免积分区间出现异常。

对高齿数齿轮，r_f可能大于r_b，那整个齿侧都是渐开线，几何处理就省心很多。实际项目里，模数3、齿数25左右的齿轮，基圆通常比齿根圆大一点，所以别默认“齿根以上全是渐开线”。

3.3 齿厚函数和积分域怎么定

在悬臂梁坐标里，把原点放在齿根处，x轴沿径向指向齿顶，载荷作用点距离齿根为d = r_load − r_f。任意截面位置x对应半径r = r_f + x。这个简化把齿廓展开当成沿径向变化，对常规直齿轮误差很小。求每个截面的半厚度h(x)，就可以算惯性矩和面积。

积分上界是载荷点半径r_load，不是齿顶，也不是整个齿高。接触点在啮合线移动时，r_load从一端的齿根区域扫向另一端的齿顶区域，所以每个啮合位置都要重新确认上界。做批量计算时，通常把离散接触点从中点附近一路布到齿顶，再反过来，组成完整啮合行程。

齿顶处齿厚最薄，截面惯性矩最小，弯曲刚度贡献也小。但要注意，接触点在齿顶附近时，齿根截面变形累积最大，所以曲线形状是两头低中间高，整个变化趋势可以从梁的长度和截面变化两个角度解释。

4. 数值实现：切片积分法

4.1 计算流程总览

整个数值流程分五步。第一步，输入齿轮副基本参数，算出各个特征圆半径。第二步，以某个齿轮转角为相位，确定当前接触点半径，或者更严格地确定沿啮合线的位置。第三步，把轮齿沿径向切成几百个薄片，逐片计算截面半厚度、惯性矩、面积。第四步，用矩形或梯形积分累加三项柔度，再加上赫兹接触柔度，得到单齿对啮合刚度。第五步，按重合度判断当前有几个齿对同时参与啮合，把所有啮合齿对的刚度并联叠加，得到整个啮合周期的TVMS曲线。

切片数一般取500到1000就够。我测试过，200片以下曲线的毛刺比较明显，500片之后积分值几乎不再变化，1000片求解时间也完全可以接受。真正影响精度的是齿根几何和角度处理，不是切片个数。

4.2 切片积分的核心代码实现

下面用Python写一个可运行的示意版本，几何参数都按毫米和兆帕输入，算出来的刚度单位是N/mm。这样直接把单位代入公式，避免换算错误。

import numpy as np # 齿轮基本参数 m = 3.0 # 模数 mm z = 25 # 齿数 alp = np.deg2rad(20.0) b = 20.0 # 齿宽 mm E = 2.06e5 # 弹性模量 MPa nu = 0.3 G = E / (2.0 * (1.0 + nu)) ha_ = 1.0 # 齿顶高系数 c_ = 0.25 # 顶隙系数 r_p = m * z / 2.0 r_b = r_p * np.cos(alp) r_a = r_p + ha_ * m r_f = r_p - (ha_ + c_) * m s_p = np.pi * m / 2.0 def involute(theta): return np.tan(theta) - theta def half_thickness(r): # r小于基圆时按基圆处齿厚延拓，避免arccos越界 r_eff = max(r, r_b * 1.0001) alpha_r = np.arccos(r_b / r_eff) s_r = r_eff * (s_p / r_p + involute(alp) - involute(alpha_r)) return s_r / 2.0 def tooth_flexibility(r_load, N=500): # 从齿根到载荷点离散 r_grid = np.linspace(r_f, r_load, N) dr = r_grid[1] - r_grid[0] d = r_load - r_f sum_b = 0.0 sum_s = 0.0 sum_a = 0.0 for r in r_grid: h = half_thickness(r) if h <= 0: continue Ix = 2.0 * b * h**3 / 3.0 Ax = 2.0 * b * h x = r - r_f # 简化取分度圆压力角；严格计算应逐点求载荷线夹角 a1 = alp sum_b += np.cos(a1)**2 * (d - x)**2 / (E * Ix) * dr sum_s += 1.2 * np.cos(a1)**2 / (G * Ax) * dr sum_a += np.sin(a1)**2 / (E * Ax) * dr return sum_b + sum_s + sum_a def mesh_pair_stiffness(r_load): flex_h = 4.0 * (1.0 - nu**2) / (np.pi * E * b) flex_tooth = tooth_flexibility(r_load) return 1.0 / (flex_h + flex_tooth)

这段代码只算了单个齿轮的单侧轮齿柔度。实际双齿轮副还要把主动轮和从动轮的几何分别算一遍，再把两个轮的柔度加在一起。从动轮齿数不同，对应半径和载荷点位置不同，函数需要重构成齿轮类结构，分别传入模数、齿数、齿宽等参数。

这里有一个工程简化：half_thickness函数在r小于基圆时按基圆处齿厚延拓，相当于把过渡区当直线近似。我项目里对比过，这样做出来的根弯曲刚度误差在可接受范围内，而且代码很稳。要是追求更精确，可以用实测或刀具参数生成齿根过渡曲线，但绝大多数TVMS分析不需要走到那一步。

4.3 单齿对啮合刚度的合成逻辑

轮齿啮合时，主动轮齿面和从动轮齿面沿法线方向接触。两个轮齿各自的弯曲、剪切、轴向压缩变形串联，赫兹接触变形也串在同一个载荷路径上，所以单对齿的总柔度写成：

1/k_pair = 1/k_h + (1/k_b1 + 1/k_s1 + 1/k_a1) + (1/k_b2 + 1/k_s2 + 1/k_a2)

计算时，接触点沿啮合线的位置会同时反映在主动轮齿和从动轮齿上。主动轮接触点半径和从动轮接触点半径不一样，需要分别求解。通常先把啮合线总长度按基节离散，得到若干接触位置，再对每个位置算两个轮的载荷点半径。

如果两个齿轮材料、模数相同，齿数不同，主动轮和从动轮的齿厚函数也不同。不要想当然地认为两轮刚度对称，尤其是齿数差别大的时候，主动轮接触点靠近齿根时，从动轮接触点可能正好在齿顶，两边柔度差异很大。

4.4 重合度与多齿对交替叠加

直齿轮副的重合度ε一般介于1到2之间。这意味着大部分时间有两对齿同时啮合，只有一小段区域是单齿对啮合。总TVMS曲线是由两个齿对刚度并联叠加出来的。比如当前时刻第一对齿的接触点在某个位置，刚度为k_A，同时第二对齿也在啮合，刚度为k_B，那整体刚度k_total = k_A + k_B。

重合度用齿轮几何直接算：

ε = [√(r_a1² − r_b1²) + √(r_a2² − r_b2²) − a sin α_p] / (π m cos α_p)

其中a = (z₁ + z₂) m / 2为中心距。这个式子的几何含义是：啮合线有效长度除以基节，就是平均同时啮合的齿对数。

实现多齿对叠加时，把驱动轮转角作为相位，按照齿距把两对齿的接触位置错开一个基节的整数倍。相位关系没对齐的话，曲线会在双齿区和单齿区交界处出现明显跳变。先画单齿对刚度曲线，再按相位错位叠加成总曲线，是最不容易出错的方式。

5. 算例：健康齿轮TVMS曲线长什么样

5.1 算例参数与基本结果

我拿一组常规参数做验证：小齿轮z₁=25，大齿轮z₂=31，模数m=3，齿宽b=20 mm，压力角20°，材料为合金钢，E=2.06×10⁵ MPa，ν=0.3。中心距a = (25+31)×3/2 = 84 mm。重合度按上式算下来在1.67附近，说明单齿区占一个基节的四成左右，其余都是双齿区。

用上面代码把啮合线离散成约100个位置，每个位置分别计算两个轮齿的柔度。单齿对刚度结果大致在1.5×10⁸到2.3×10⁸ N/m之间波动，折算到整个齿宽上就是每毫米齿宽大约7500到11500 N/mm。这个量级和文献中同参数齿轮的解析结果对得上，说明势能法得到的结果是可靠的。

5.2 双齿区和单齿区的衔接特征

一个完整啮合周期里，总刚度曲线会呈现明显的“哑铃形”起伏。双齿区因为有两条载荷路径，总刚度偏高；切入单齿区后，只剩一对齿承载，刚度陡然下降，形成局部谷值。进入下一组双齿区时，刚度又抬升。由于齿轮不停旋转，这种高低交替就形成周期性波动，波动频率就是轮齿啮合频率。

更仔细看单齿对刚度曲线本身，接触点刚进入啮合时往往靠近某一齿的齿根，悬臂长度短、根部截面厚，刚度大；随着啮合点往齿顶移动，悬臂变长、齿厚变薄，刚度下降；到脱离前又会快速回升。把这条单齿对曲线与多齿对叠加逻辑结合，总TVMS曲线的凹凸位置就能解释清楚。

实际我在曲线里观察到的最小值点并不是恰好位于单齿区正中，而是偏向单齿区与双齿区交界处。原因是单齿区前后两个齿对在不同位置分担载荷，重叠出来刚度的相对大小并不对称。做故障诊断时，这个最小值点的转角位置要标定准确，裂纹定位才可靠。

5.3 验证、误差来源与收敛性

验证时我会做两个检查：第一，收敛性检查，切片数从200、500到1000，看刚度曲线之间的最大偏差是否控制在1%以内；第二，交叉验证，对同样的齿轮副用有限元静力接触算几个特征相位，比较单齿对刚度偏差。常规设计下，势能法和有限元的偏差通常在5%以内，齿根过渡区处理粗糙时可能拉到10%。

误差主要来自三块：齿根过渡曲线近似、载荷角α₁固定为分度圆压力角、忽略轮体柔度。前三者对曲线中段影响小，对齿根附近影响最大。如果做的是齿根裂纹趋势分析，只要基线和带裂纹模型用同一套近似，相对变化趋势仍然可信，没必要为了绝对精度强行加复杂几何。

6. 实操避坑与常见问题

6.1 单位制不统一是最大的坑

我见过太多次结果差三个数量级的问题，最后都是单位制搞混。刚度公式里如果力用N、长度用mm，那弹性模量必须用MPa（N/mm²），惯性矩单位是mm⁴，面积单位是mm²，积分出来的柔度单位是mm/N，取倒数得到N/mm。有人习惯把弹性模量写成2.06×10¹¹ Pa，长度又用mm，最后差出10⁶倍。最稳妥的办法是全部用mm和MPa，代码里加注释，最后再统一检查一遍。

6.2 齿根圆和基圆边界处理不当

前文已经强调过r_f和r_b的关系。程序里如果直接对r < r_b做arccos，会得到NaN或者虚数；即使强行clip，也可能生成不存在的齿厚分布。碰到这种情况，第一选择是用真实齿根过渡曲线，第二选择是近似延拓，但要在结果说明里讲清楚。不要为了曲线好看把齿根截面积分乱改，否则后面拿它做故障判断会误导自己。

6.3 载荷角α₁要不要逐点算

我在初版计算里直接用了分度圆压力角，结果曲线形态正确，但峰值略偏高。后来把载荷角做成逐点几何插值后，整体刚度小幅下降，曲线过渡更平滑。做法是：根据啮合线方向、主动轮转角、接触点半径，算出法向力方向与该轮齿对称线的夹角，然后代入cos²项和sin²项。这个修正对弯曲刚度影响约10%，对轴向刚度影响更大但轴向刚度本身占比小，所以总体影响有限。

6.4 多齿对叠加的相位对不齐

多齿对叠加看似简单，实际非常容易差半个齿距。建议先把主动轮一个齿距等分成若干相位格点，将第一对齿的接触半径计算好，然后让第二对齿的接触半径相对第一对偏移一个基节对应的角度。基节是π m cos α_p，对应角度用基圆半径折算。画图时把横轴统一换算成接触线位移或转角，不要混用分度圆弧长和基圆弧长。

6.5 曲线出现尖角或负刚度

尖角通常来自切片数太少或齿厚函数在基圆处不连续。负刚度几乎一定是半厚度h(x)出现了负值或者积分方向反了。排查时把half_thickness函数单独输出，检查它在齿根到齿顶区间是否单调递减，正常情况半齿厚应该从齿根向齿顶逐渐变小。如果看到突变，先确认齿根圆和基圆的判断逻辑。

7. 这套结果后续能怎么用

算完健康齿轮TVMS，我给自己的经验是别急着收工。可以在同一套代码里继续做三件很有价值的事。第一，修形模拟：把齿廓上的接触点位置做微量修形偏移，重新计算接触半径和齿厚分布，就能对比修形前后刚度突变量。第二，裂纹模拟：在齿根截面引入一个局部开口深度，弯曲刚度积分时把惯性矩做局部折算，能很快看出TVMS在哪个转角位置塌陷，这也是不少故障诊断论文里做灵敏度分析的标准起手式。第三，动力学耦合：把计算得到的TVMS曲线做成傅里叶级数形式，代入集中质量齿轮动力学模型，可以复现边带频谱和振动响应调幅现象。

实际操作中我也踩过几次坑，最大的体会是：方法传统不等于粗糙，关键是每一步近似都要知道自己在近似什么。势能法看起来只有几个积分式，但几何建模和相位关系才是真正拉开差距的地方。算健康齿轮基线的过程，也是把整个齿轮啮合过程重新梳理一遍的过程，后面再引入故障、修形、误差，思路都会清晰很多。你一旦把这条基线跑通了，从“会算”到“能分析”的距离其实已经跨过一大半。