我做齿轮动力学分析有一阵子了，最开始就是用商业软件算箱体振动，但每次都被黑箱式的过程搞得一头雾水：输出一大堆模态参与因子，现场工程师真正想知道的“齿频附近为什么会有共振峰”“把齿宽加厚能不能把振动压下去”这类问题，软件讲不清楚。后来我干脆自己写了一套齿轮动力学求解程序，从数学模型、数值积分到后处理全部自己控制，反而把很多原理彻底弄明白了。这篇就把这套程序的原理和实战经验拆开讲：怎么建模型、怎么选求解器、哪些参数最容易出错、怎么验证程序算得对。适合刚接触齿轮动力学仿真的机械工程师、研究生，以及所有想摆脱“只会点按钮”状态的人。
1. 为什么需要自己写齿轮动力学求解程序——从工程痛点说起
1.1 设计手册算不清的问题：动态载荷与共振
齿轮传动的经典设计流程里，我们通常用齿根弯曲强度和齿面接触强度去校核静强度，公式把载荷当成一个静态值乘以工况系数。但这个做法有个很大的盲区：齿轮啮合过程中，单双齿交替、啮合刚度的周期性变化、齿轮误差的周期激励，都会让实际动态载荷明显大于静态载荷。尤其当啮合频率或其倍频与系统固有频率重合时，动载系数可以拉到2倍甚至更高，这时候静强度校核全过也没用，箱体该裂还是裂。
我遇到过一个很典型的案例（某模拟项目X）：一个二级减速器在某个转速区间出现明显的齿轮敲击声，静强度校核完全合格，换了好几家润滑油都没解决。最后溯源到输入级主动轮的转频倍频与轴系一阶扭转固有频率接近，产生扭转共振。这种问题靠查手册的公式是查不出来的，必须把“齿轮副-轴-轴承-箱体”作为一个动力学系统去算。
1.2 商业软件之外的自研路线：谁适合这么做
很多人会问：ANSYS、Simpack、Romax这类工具不是能算吗？能算，但有几个现实问题。第一是费用，正版授权对一个中小型团队来说不便宜。第二是建模过程不透明，软件内部的接触算法、阻尼假设固化了，你很难精确控制某个参数，比如啮合阻尼比。第三是最要命的——计算速度。完全有限元模型算一个工况的瞬态响应可能要跑十几个小时，参数优选根本没法做。
自研求解程序的价值就在这：把动力学模型降成集中参数（集中质量、弹簧、阻尼），单工况几秒钟就算完，批量扫参数特别方便。虽然精度上限不如精细有限元，但在工程前期评估和趋势预测这个层面，性价比极高。这里的核心词是“齿轮动力学求解程序”，它本质上是“用尽量简单但物理正确的模型，把齿轮系统的振动响应算出来”的一套代码。适合的人群，是那些已经具备基础振动理论、同时又愿意花时间把物理模型吃透的工程师和研究者。
2. 动力学模型的建立：从齿轮副到多自由度系统
2.1 单自由度啮合线模型：用等效质量把问题压简
齿轮动力学里最经典的入门模型，是把一对齿轮副化简成沿啮合线方向的单自由度振动系统。思路是这样的：齿轮啮合时，作用力沿啮合线传递，齿面接触变形也发生在啮合线方向，所以可以把两个旋转自由度映射到啮合线位移上。
定义啮合线位移为
[ x = r_{b1}\theta_1 - r_{b2}\theta_2 ]
其中 (r_{b1})、(r_{b2}) 是两个齿轮的基圆半径，(\theta_1)、(\theta_2) 是角位移。这样定义的好处是：当 (x>0) 表示齿面脱开方向的相对位移，(x<0) 表示齿面压紧。系统的等效质量由两个齿轮的转动惯量折算到啮合线方向得到，注意这里的正确表达式是
[ m_e = \frac{I_1}{r_{b1}^2} + \frac{I_2}{r_{b2}^2} ]
我拿一个常用参数举例：主动轮模数2mm、齿数23、基圆半径约21.6mm，转动惯量取 (2.25 \times 10^{-4}\ \mathrm{kg\cdot m^2})，从动轮齿数46、基圆半径约43.2mm，转动惯量约 (9 \times 10^{-4}\ \mathrm{kg\cdot m^2})，算出来等效质量大约1kg。这个量级很有用——你可以立刻估算固有频率：平均啮合刚度 (k_0) 如果是 (2.5 \times 10^8\ \mathrm{N/m})，固有频率就是 (f_n \approx \frac{1}{2\pi}\sqrt{k_0/m_e} \approx 2500\ \mathrm{Hz})。看到这个数，你就能快速判断工作转速下有没有共振风险，整个过程只需要按几下计算器。
2.2 刚度激励与误差激励：两个最核心的周期源
齿轮系统为什么会振动？两个最核心的周期激励源：刚度激励和误差激励。
刚度激励来自啮合过程中单双齿交替。重合度通常在1到2之间，意味着有时候一对齿啮合，有时候两对齿同时啮合，综合啮合刚度跟着周期变化。这个变化的基频就是啮合频率：
[ f_z = \frac{n \cdot z}{60} ]
其中 (n) 是转速（r/min），(z) 是齿数。比如3000r/min、23齿，啮合频率就是1150Hz。时变啮合刚度写成：
[ k_m(t) = k_0 + \sum_{i=1}^{\infty} k_i \cos(i \omega_z t + \varphi_i) ]
实现时一般取前3到5阶谐波就够了，因为高阶分量幅值迅速衰减。
误差激励来自齿轮的制造误差、安装误差、齿距偏差和齿形偏差。这些误差导致齿轮在匀速啮合时实际传动比周期性波动，等效成一个位移激励 (e(t))，作用在啮合线方向上。最终的运动方程写成：
[ m_e \ddot{x} + c_m \dot{x} + k_m(t) x = F_d + k_m(t) e(t) ]
右边第一项是外载荷，第二项是误差激励。这里要格外注意：误差项是乘上啮合刚度的，因此齿越“硬”，同样的误差激起的动态力越大。
2.3 多自由度扩展：轴承支撑、轴系扭振与箱体耦合
单自由度模型适合理解机理和快速估算，但工程上不够用。轴的弯曲变形、轴承支撑刚度不对称、两个齿轮轴的扭转耦合，都会显著影响响应。我实际用的程序是一个多自由度集中参数模型，广义坐标包括：
- 每个齿轮的扭转角位移
- 每根轴在轴承位置的横向位移
- 必要时加上齿轮的摆动角
系统方程为：
[ \mathbf{M}\ddot{\mathbf{q}} + \mathbf{C}\dot{\mathbf{q}} + \left[\mathbf{K}_b + \mathbf{K}_m(t)\right]\mathbf{q} = \mathbf{F}(t) ]
其中 (\mathbf{K}_b) 是轴承和轴的支撑刚度矩阵，(\mathbf{K}_m(t)) 是时变啮合刚度矩阵。啮合刚度在小位移假设下可以线性化，但当齿面脱啮或发生齿背接触时，刚度会变成非线性的，程序里需要加判断分支。多自由度模型的优势在于能算出每个测点的振动响应，直接和试验对标；代价是参数辨识难度上升，轴承刚度、轴段刚度都得想办法估准。
3. 数值求解方法选型：直接积分、模态叠加与谐波平衡
3.1 四阶Runge-Kutta：原理和步长控制
动力学方程本质上是一个二阶常微分方程组。最直观的做法是化成状态空间形式，用四阶Runge-Kutta（RK4）做时间推进。RK4的迭代公式如下：
[ \begin{aligned} k_1 &= f(t_n, y_n) \ k_2 &= f(t_n + \frac{h}{2}, y_n + \frac{h}{2} k_1) \ k_3 &= f(t_n + \frac{h}{2}, y_n + \frac{h}{2} k_2) \ k_4 &= f(t_n + h, y_n + h k_3) \ y_{n+1} &= y_n + \frac{h}{6}(k_1 + 2k_2 + 2k_3 + k_4) \end{aligned} ]
RK4实现简单，数值耗散也小，适合中低自由度系统。但它是显式方法，步长受稳定性和精度限制。经验规则是：系统最高固有频率对应的周期至少要采20个点以上。比如前面算出的固有频率2500Hz，步长就要小于 (1/(2500 \times 20) = 20\ \mu s)。如果扫频到高速段，啮合频率高、刚度又时变，步长还得进一步缩小。
所以我一般用自适应步长：每步估算局部误差，误差超出容差就自动减半步长，误差过小就翻倍步长。这样可以避免人为拍脑袋定步长，也能处理刚度激励突变处的瞬态。
3.2 Newmark-β：为什么工程上更常用
如果系统自由度多起来，比如有十几个坐标，再叠加齿轮啮合的强非线性，纯显式RK4的时间步长会被压得很小，计算时间太难看。这时候我推荐用Newmark-β法。它的核心思想是通过加速度线性外推把位移和速度递推出来，参数β控制加速度在时间步内的加权方式：
[ \begin{aligned} \dot{q}_{n+1} &= \dot{q}n + (1-\gamma) h \ddot{q}n + \gamma h \ddot{q}{n+1} \ q{n+1} &= q_n + h \dot{q}_n + \left(\frac{1}{2}-\beta\right) h^2 \ddot{q}n + \beta h^2 \ddot{q}{n+1} \end{aligned} ]
取平均加速度法的参数组合 (\gamma = 0.5)、(\beta = 0.25)，该方法是无条件稳定的，步长可以放宽到不再受最高固有频率的限制，只需要保证能捕捉激励的高频成分即可。这在工程仿真里特别香——齿轮啮合频率往往是几百到几千赫兹，但系统刚度时变带来的非线性成分可能到更高频，Newmark-β能用一个相对稳定的步长覆盖整个频段。
但无条件稳定不等于无条件准确，高频模态的响应会被平均掉，造成相位失真。所以我在共振分析这种关注精确幅值和相位的场景，会回到RK4或者更精细的步长设置。
3.3 什么时候才能偷懒用频域方法
齿轮系统有个特点：激励源几乎都是周期性信号，理论上稳态响应也是周期的。因此可以用谐波平衡法或者傅里叶级数展开，把时域方程变成频域代数方程来解。好处是避免长时间积分，算一个工况可能只要零头时间。
条件限制也明显：系统必须是弱非线性的。齿轮脱啮、齿背接触、轴承间隙这些强非线性一旦发生，频率响应里会出现连续频谱，谐波平衡法就力不从心了。我的做法是：先跑一次快速时域仿真判断啮合状态，如果始终是“单齿-双齿”交替且无脱啮，可以用谐波平衡法提速；一旦监测到位移超过齿侧间隙，立刻转回时域直接积分。模块化程序的好处就在这——求解器可以按工况切换。
4. 程序实现的关键模块：架构、数据流与一段可跑示例
4.1 模块划分与数据结构：避免全局变量地狱
程序写大了以后最容易乱的是参数传递。我的经验是分四个模块，每个模块只干一件事：
params：集中存放所有几何参数、材料参数、工况参数，统一用SI单位mesh：计算时变啮合刚度、传递误差激励，输出随时间变化的激励项dynamics：组装质量、阻尼、刚度矩阵，定义状态方程右端函数postprocess：时域响应、FFT频谱、RMS值、瀑布图生成
数据结构上用dataclass或者嵌套字典保存参数，避免到处写全局变量。例如齿轮副参数可以定义成：
from dataclasses import dataclass @dataclass class GearPair: z1: int # 主动轮齿数 z2: int # 从动轮齿数 mn: float # 法向模数，单位m alpha_n: float # 法向压力角，单位rad r_b1: float # 主动轮基圆半径 r_b2: float # 从动轮基圆半径 I1: float # 主动轮转动惯量 I2: float # 从动轮转动惯量这样后面所有模块只需要传入这个对象，不会出现“函数里改了某个变量，另一个函数跟着出错”的连锁反应。
4.2 最小可运行示例：单自由度RK4求解
我给出一个能直接跑通的单自由度模型代码，便于你建立第一版程序。这段代码包含刚度激励和误差激励，采样时长取到0.4s，确保系统进入稳态：
import numpy as np class SDOFGear: def __init__(self, m_e, k0, k1, z, rpm, zeta, F0): self.m_e = m_e self.k0 = k0 self.k1 = k1 self.z = z self.omega_z = 2.0 * np.pi * z * rpm / 60.0 self.zeta = zeta self.F0 = F0 self.c = 2.0 * zeta * np.sqrt(k0 * m_e) def stiffness(self, t): # 时变刚度，取基频分量近似 return self.k0 + self.k1 * np.cos(self.omega_z * t) def error(self, t): # 传递误差简化模型，幅值5微米的简谐波 return 5e-6 * np.sin(self.omega_z * t) def rhs(self, t, y): x, v = y k = self.stiffness(t) e = self.error(t) a = (self.F0 + k * e - self.c * v - k * x) / self.m_e return np.array([v, a]) def rk4_step(f, t, y, h): k1 = f(t, y) k2 = f(t + 0.5*h, y + 0.5*h*k1) k3 = f(t + 0.5*h, y + 0.5*h*k2) k4 = f(t + h, y + h*k3) return y + (h/6.0)*(k1 + 2.0*k2 + 2.0*k3 + k4) # 参数：等效质量1kg，平均刚度2.5e8 N/m，刚度波动0.15倍 gear = SDOFGear(m_e=1.0, k0=2.5e8, k1=3.8e7, z=23, rpm=3000, zeta=0.05, F0=1200.0) t0, t_end, dt = 0.0, 0.4, 1e-5 steps = int((t_end - t0) / dt) t_arr = np.empty(steps) x_arr = np.empty(steps) y = np.array([0.0, 0.0]) t = t0 for i in range(steps): t_arr[i] = t x_arr[i] = y[0] y = rk4_step(gear.rhs, t, y, dt) t += dt注意步长我选了10微秒，对应最高2500Hz固有频率一个周期40个点，精度足够。跑完可以观察稳态段的响应幅值，对比静态变形 (F_0/k_0 = 1200/2.5e8 = 4.8\ \mu m)，动态幅值会比这个大，大多少就取决于阻尼和激励频率的接近程度。把这套流程跑通，再扩展多自由度就只是矩阵维度问题。
4.3 后处理：时域、FFT与瀑布图的正确打开方式
算完响应不算完，工程判断靠后处理。我在程序里最常用的三种输出：
时域曲线主要看瞬态衰减和稳态幅值，判断是否发生脱啮——如果啮合线位移 (x) 大于齿侧间隙，说明齿面脱开了，程序输出里会直接显示标志位。
FFT频谱看频率成分。用np.fft.rfft对稳态段加汉宁窗再变换，横轴标成频率，重点关注啮合频率及其倍频、边带（调制成分往往是故障特征）。注意采样时长要覆盖至少10个啮合周期，否则频域分辨率不够，边带根本分不出来。
瀑布图用于扫频工况：横轴是频率、纵轴是转速、颜色或高度是幅值。一张图就能看出当转速变化时共振峰怎么移动、和哪阶固有频率交汇。画瀑布图的核心是必须保证每个转速点都进入稳态后再采样，否则瞬态分量会污染频谱。
5. 实测中最容易踩的三个坑：阻尼、符号、初始冲击
5.1 阻尼参数混乱：临界阻尼、相对阻尼和瑞利阻尼
阻尼是最难估的参数，也是仿真和试验差别的头号来源。做齿轮动力学时，最容易踩的坑是把临界阻尼系数和阻尼比搞混。临界阻尼系数 (c_c = 2\sqrt{km})，而运动方程里和速度相乘的 (c_m) 通常是 (c_m = 2\zeta\sqrt{k m_e})，其中 (\zeta) 是阻尼比。齿轮啮合副的阻尼比经验值大概在0.03到0.10之间，千万不能用干摩擦或者结构阻尼的经验值去套。
多自由度系统用瑞利阻尼 (C = \alpha M + \beta K)。取两个关心的振型频率 (\omega_1)、(\omega_2) 和对应的阻尼比 (\zeta_1)、(\zeta_2)，解二元一次方程：
[ \begin{aligned} \alpha + \beta \omega_1^2 &= 2 \zeta_1 \omega_1 \ \alpha + \beta \omega_2^2 &= 2 \zeta_2 \omega_2 \end{aligned} ]
这里有个隐性约束：(\omega_1) 和 (\omega_2) 不能倍频太近，否则矩阵会病态。我在实际程序里用这个公式前会先画一下阻尼比-频率曲线，确认感兴趣的频段内阻尼比平滑合理。
5.2 传递误差符号与激励相位
误差激励项的符号是我调试程序时踩过最深的坑。方程式里误差项是 (k_m(t) e(t))，如果符号反了，激励相位差180度，共振峰的位置、幅值全部对不上。很多刚接触的人以为反正激励是周期的，符号无所谓，实际影响很大。
怎么检查？把传递误差设成单频简谐波，激励频率固定在某阶固有频率附近，算稳态响应。理论规律是：低于共振点时，位移响应与激励力同相；高于共振点时，响应滞后180度。程序输出做了相位谱，直接比对这个规律，符号错了立刻暴露。任何程序升级后，我都会跑一遍这个相位检查，三分钟内就能确认建模没破坏。
5.3 启动冲击与初始条件的静力学预处理
另一个新手必踩的坑：给恒定扭矩直接开始积分。这等于对系统施加了一个阶跃载荷，响应里会叠加严重的瞬态冲击，频谱上出现宽频隆起，严重影响共振峰的判断。
我常用的解决办法有两种。一种是对外力做斜坡处理，比如前0.1秒让外载荷从0线性升到目标值：
[ F_d(t) = F_0 \cdot \min(1, t / T_{\rm ramp}) ]
另一种更正规：先求解静平衡，把静变形作为初始条件。用 (K(0) q_0 = F) 解出初始位移，初始速度设为0。这样启动那一瞬间系统已经处于平衡位置，载荷冲击基本消除。两条路可以一起用——斜坡主要抑制载荷冲击，静预处理主要抑制刚度突变。
6. 结果验证与工程落地：怎么让人相信你的程序算得对
6.1 退化测试：把时变刚度设成常数
程序写完了先别急着算复杂工况，第一步是退化测试。把时变刚度 (k_m(t)) 设成常数，误差激励设成零，模型退化成一个标准线性单自由度系统。此时如果外载荷是恒定的，系统最终应停在静变形位置，动态响应衰减归零。我常用来验证：算0.2秒，看位移是否收敛到 (F_0/k_0)，误差在1%以内说明积分器、矩阵组装没有低级问题。
如果退化模型还带简谐激励，那就更严格了。线性系统有解析解：稳态振幅为
[ |X| = \frac{F_0/k}{\sqrt{(1-\beta^2)^2 + (2\zeta\beta)^2}} ]
其中 (\beta = \omega / \omega_n)。用程序算出来的稳态幅值和这个解析解对比，幅值差应在0.5%以内（若步长足够），相位差也应在允许范围内。这一步过了，说明核心积分过程是可信的。
6.2 与线性解析解对拍：一个共振点一个共振点地验
退化测试之后，我习惯扫一版频率响应。在恒定激励幅值下，从0.5倍固有频率扫到2倍固有频率，把每个频率点的稳态响应幅值画出来。和解析解曲线重在图里叠加，你会看到峰值基本叠合，这时心里才踏实——共振频率、幅值、带宽都对了。
这个过程中有个细节：扫频时每个频率点都要跳过瞬态段再采样。我一般取30个激励周期作为稳定时间，再取后续10个周期做幅值平均。否则扫出来的峰值会偏低，因为采样时刻响应还没真正进入稳态。这个细节是纯从实操里学到的，文档里很少写。
6.3 台架实测对标：误差能压到多少才算可用
最后一步是拿真实齿轮箱的测试数据对标。我在某模拟项目X的减速器上做过一次验证：输入轴轴承座贴三向加速度传感器，转速从1200r/min升到3600r/min，记录箱体振动；同时用程序预测轴承座附近的振动响应。实测频谱里，啮合频率及其边带的分布位置、共振转速区间，和程序预测基本一致，共振峰对应的转速点误差在5%以内；但齿轮啮合基频的绝对幅值，程序预测值和实测差了大概18%。
这个偏差主因是阻尼比标定不足、轴承等效刚度简化。要压到10%以内，就得用锤击模态试验反辨识阻尼，并且把轴承刚度的非线性特性标出来。但即便有18%的幅值误差，程序做参数趋势研究是完全够用的——比如齿宽增加10%振动能降多少、改变支承位置共振转速往哪移，这些相对变化量程序算得很稳。工程上很少需要绝对幅值有多准，大多数决策都是相对比较。这也是为什么我一直推荐用自研程序做前期选型和故障复盘。
我个人做这套程序的体会是：先跑通单自由度模型，别一上来就堆矩阵自由度和商业级后处理。单自由度模型能帮你把刚度激励、误差激励、阻尼、步长这些基本概念摸清楚，后面加自由度、加非线性分支都是顺理成章的事。验证一步也不能省，每加一个模块就退化测试一次。等程序算出来的趋势在台架上被验证过一轮，你就有底气用它指导实际产品设计——这种感觉，是黑箱软件给不了的。
