齿轮啮合刚度傅立叶级数展开：从频谱边频带到谐波参数提取

齿轮箱振动超标，频谱图上除了啮合频率及其倍频，两侧还冒出一堆边频带——这几乎是每一个做齿轮传动系统设计的人都撞过的墙。我当时排查了很久，最后才把问题锁定在时变啮合刚度上。说白了，齿轮在啮合过程中，参与啮合的轮齿对数在周期性变化，系统刚度也跟着周期性波动，这就是齿轮振动噪声最核心的激励源。而要定量分析这个激励，把啮合刚度展开成傅立叶级数，几乎是绕不开的一步。

这篇文章我就围绕“齿轮啮合刚度傅立叶级数展开程序”这件事，完整梳理一下背后的物理逻辑、程序实现思路、关键代码段，以及我实际踩过的一些坑。无论你是做齿轮箱NVH的工程师、研究传动系统动力学的学生，还是刚接触齿轮程序化分析的入门者，这篇文章都能帮你把这块硬骨头啃下来。

1. 齿轮啮合刚度的时变特性：为什么它天生就是周期函数

1.1 刚度的物理来源：单双齿交替啮合的“呼吸效应”

齿轮啮合刚度，通俗讲是轮齿抵抗弹性变形的能力，单位常用N/(mm·μm)或者N/mm。它不是一个恒定值，而是随啮合位置变化的。原因很简单：齿轮传动时，啮合过程是单齿啮合区和双齿啮合区交替出现的。

在双齿啮合区，两对轮齿同时分担载荷，系统表现出的综合刚度高；在单齿啮合区，只有一对轮齿承载，刚度明显下降。这个高低交替在时域上形成了类似方波或梯形的周期波动，我习惯管它叫“呼吸效应”——因为刚度曲线看起来就像在一呼一吸地起伏。

重合度是决定这个波形长相的关键参数。重合度在1到2之间时，一个啮合周期内会有一段单齿区加一段双齿区；重合度大于2时，会出现三齿啮合区，波形形态又会不一样。设计齿轮时，重合度越高，刚度波动幅度相对越小，传动越平稳，这是后话。

1.2 为什么时变刚度是振动噪声的“元凶”

从动力学角度看，齿轮副可以简化成一个参数激励系统：

m·x'' + c·x' + k(t)·x = F(t)

这里k(t)就是时变啮合刚度。这个方程和普通强迫振动最大的区别在于：刚度项本身就是时间的函数。哪怕外载荷F(t)恒定不变，系统也会因为刚度的周期性变化而被“激励”起来，产生参数振动。

所以你会看到，齿轮箱振动频谱上，啮合频率fm及其整数倍频处总有明显峰值，即使负载非常平稳。fm = z·n/60，z是齿数，n是转速（r/min）。啮合频率的本质，就是轮齿啮合事件重复发生的频率，而刚度波形每个周期重复一次，它的基频正好等于fm。

既然刚度是周期函数，那就能用傅立叶级数展开。展开之后，时变刚度就变成：

k(t) = k₀ + Σ|K_m|·cos(2π·m·fm·t + φ_m)

k₀是平均刚度，决定系统的固有频率；后面那些谐波项才是真正的激励源，每一项的幅值|K_m|和相位φ_m，直接影响齿轮箱在对应频率上的振动响应。这就是为什么“齿轮啮合刚度傅立叶级数展开程序”在工程里那么重要——它把一条复杂的刚度曲线，压缩成一组简洁的谐波参数，方便后续做任何动力学分析。

2. 程序实现的核心思路：从刚度曲线到谐波参数

2.1 刚度曲线的获取途径

写程序第一步是拿到一个啮合周期内的刚度数据。实际工程中有三种主流途径：

• 解析公式法：基于石川公式、Weber-Banashek公式或者势能法，根据齿轮几何参数直接计算不同啮合位置的刚度值。速度快、参数化方便，适合方案设计阶段批量扫描。缺点是精度一般，对修形、微观几何的考虑比较粗糙。
• 有限元法：建立轮齿接触模型，通过逐步旋转齿轮副获取不同啮合位置的接触刚度和弯曲刚度。精度高，能考虑齿廓修形、齿向修形、甚至实际接触斑点，但计算量感人，一个完整啮合周期的扫描动辄几十个小时。
• 试验测试法：用应变片或激光测振仪实测轮齿变形。最真实，但成本高、周期长，通常只用于验证阶段。

从程序设计角度讲，无论数据来自哪种途径，后面傅立叶展开部分的处理逻辑是完全一致的。我一般建议前期用解析法或者简化模型先把程序链路跑通，确认傅立叶展开模块的精度和稳定性，再接入有限元结果，否则一上来就用有限元数据调试程序，排查问题会非常痛苦。

2.2 离散化与序列预处理

假设我们已经拿到了一个完整啮合周期内的刚度数据，N个等间距采样点，记为k[0], k[1], ..., k[N-1]。在送入傅立叶变换之前，一定要做两件事：

第一，确认序列确实是整周期截断的。也就是说，k[0]和k[N-1]应该分别对应啮合周期起点和终点前的最后一个点，k[N]理论上应该等于k[0]。如果截断位置不对，或者点数没掐准，后面拉出来的频谱全是泄漏，谐波系数会严重失真。这个地方我吃过亏，后面专门讲。

第二，理解序列的直流分量。对刚度序列做FFT后，第0个bin对应的是平均刚度k₀，也就是直流分量。很多初学者上来就把FFT结果直接当激励幅值用，忽略了直流分量才是决定系统固有频率的关键参数，谐波项只负责“扰动”。程序里最好把直流和波动分开存储，后面各有各的用场。

2.3 傅立叶系数提取的数学原理

对离散序列做DFT的标准公式：

K_m = Σ_{n=0}^{N-1} k[n]·e^{-j·2π·mn/N}，m = 0, 1, ..., N-1

K_m就是第m阶谐波的复系数。实际程序里直接用FFT算法，复杂度从O(N²)降到O(N·logN)，几百上千个点都是毫秒级的事。

这里有一个关键的物理对应关系需要讲清楚。刚度曲线的基频是啮合频率fm，所以：

• K₀对应0Hz，即平均刚度；
• K₁对应频率fm，即第一阶啮合谐波；
• K₂对应频率2fm，即第二阶啮合谐波；
• K_m对应频率m·fm，即第m阶啮合谐波。

FFT输出的每个bin，对应的实际频率是：

f_bin = bin_index / N · f_s

其中f_s是采样频率。但因为刚度曲线的周期是一个啮合周期T_m = 1/fm，我们在一整个周期内采样N个点，所以等效的采样频率就是f_s = N·fm。代入前面的公式你会发现，bin_index恰好等于谐波阶次m时，f_bin = m·fm。这就是为什么“整周期采样 + FFT”时，频谱里的bin序号可以直接当作谐波阶次来读，非常直观。

2.4 程序主流程设计

整个程序的主流程我建议这样组织：

1. 输入齿轮参数（齿数、模数、压力角、变位系数、重合度等），选定刚度模型；
2. 计算一个啮合周期内的离散刚度序列k[n]；
3. 对序列做预处理（去趋势、整周期确认）；
4. 调用FFT得到复数谱；
5. 提取各阶谐波的幅值和相位；
6. 用前M阶谐波重建刚度曲线，与原始曲线对比验证；
7. 输出平均刚度、各阶谐波幅值相位表，供下游动力学分析使用。

这个流程很朴素，但每一步都有坑。接下来我直接上代码和实际算例，一步步拆开讲。

3. 傅立叶展开程序的代码实现：核心函数与算例验证

3.1 构造测试用的刚度序列

为了聚焦傅立叶展开这个主题，我这里用一个典型的单双齿交替刚度模型来生成测试序列。实际工程中，你可以把这一整段替换成自己的有限元结果或实测数据。

import numpy as np import matplotlib.pyplot as plt # 齿轮基本参数 z1 = 26 # 主动轮齿数 z2 = 58 # 从动轮齿数 m = 3.0 # 模数 mm alpha = 20.0 # 压力角 deg epsilon = 1.65 # 重合度 # 啮合刚度典型值（经验范围，单位 N/mm/μm） k_single = 14.0 # 单齿啮合区刚度 k_double = 22.0 # 双齿啮合区刚度 # 一个啮合周期内采样点数 N = 1024 # 构造刚度波形：双齿区占比由重合度决定 # 重合度1.65意味着65%的时间处于双齿区，35%处于单齿区 double_ratio = epsilon - 1.0 single_ratio = 2.0 - epsilon n_double = int(N * double_ratio) n_single = N - n_double # 为了模拟真实刚度曲线的光滑过渡，用分段余弦过渡代替理想方波 t = np.linspace(0, 1, N, endpoint=False) # 理想梯形波：上升沿/下降沿各占5%周期 transition = int(0.05 * N) k_series = np.zeros(N) for i in range(N): if i < transition: # 上升沿：从单齿刚度到双齿刚度 k_series[i] = k_single + (k_double - k_single) * (0.5 - 0.5 * np.cos(np.pi * i / transition)) elif i < transition + n_double - transition * 2: # 双齿平台区 k_series[i] = k_double elif i < n_double: # 下降沿 idx = i - (n_double - transition) k_series[i] = k_double - (k_double - k_single) * (0.5 - 0.5 * np.cos(np.pi * idx / transition)) else: # 单齿区 k_series[i] = k_single

代码里用余弦过渡替代理想方波，是为了模拟真实轮齿啮合过程中的弹性变形过渡。实际轮齿进入啮合和退出啮合时，刚度变化不是瞬时的，而是有一个渐变过程。过渡区取周期的5%左右，已经是一个比较合理的近似。

3.2 FFT谐波提取与幅值相位还原

接下来是核心的傅立叶展开部分：

def fourier_expansion(k_series, N_harmonics): """ 对刚度序列做傅立叶展开，提取前N_harmonics阶谐波 返回：直流分量k0，各阶幅值数组，各阶相位数组 """ N = len(k_series) # 做FFT K = np.fft.fft(k_series) # 归一化：DFT的幅值需要除以N才是真实的波形幅值 K_normalized = K / N # 直流分量 k0 = np.real(K_normalized[0]) # 提取前N_harmonics阶谐波幅值和相位 amplitudes = np.zeros(N_harmonics) phases = np.zeros(N_harmonics) for m in range(1, N_harmonics + 1): amplitudes[m-1] = 2.0 * np.abs(K_normalized[m]) phases[m-1] = np.angle(K_normalized[m]) return k0, amplitudes, phases # 提取前20阶谐波 k0, amps, phis = fourier_expansion(k_series, 20) print(f"平均刚度 k0 = {k0:.4f} N/mm/μm") print("前10阶谐波幅值：") for i, a in enumerate(amps[:10]): print(f" 第{i+1}阶: {a:.4f}")

这段代码有几个细节值得注意。

第一，FFT归一化。Python的np.fft.fft直接输出的是复数谱，数值大小和采样点数N成正比。要还原真实幅值，必须除以N。这是最容易犯的错误，我见过不少同事在这里栽过。

第二，单边谱的处理。实数序列的FFT结果是对称的，正频率和负频率各占一半。提取幅值时，需要对正频率部分乘以2（直流分量除外），才能得到真实的余弦波幅值。上面的代码里amplitudes[m-1] = 2.0 * np.abs(K_normalized[m])，这个2.0就是干这个用的。

3.3 用谐波重建刚度曲线验证精度

傅立叶级数展开完不完整，不能光看系数好看，必须要重建回时域曲线，和原始波形对比，才能发现问题。重建代码如下：

def reconstruct_series(k0, amps, phis, N): """ 用直流分量和前M阶谐波重建刚度序列 """ n = np.arange(N) k_recon = np.full(N, k0) for i, (a, p) in enumerate(zip(amps, phis)): m = i + 1 k_recon += a * np.cos(2 * np.pi * m * n / N + p) return k_recon # 使用前10阶谐波重建，对比原始波形 k_recon = reconstruct_series(k0, amps[:10], phis[:10], N) # 计算均方根误差 rms_error = np.sqrt(np.mean((k_series - k_recon)**2)) rms_value = np.sqrt(np.mean(k_series**2)) relative_error = rms_error / rms_value * 100 print(f"前10阶重建的相对误差: {relative_error:.4f}%") # 用不同谐波阶次验证收敛性 for M in [1, 3, 5, 10, 20]: k_recon_M = reconstruct_series(k0, amps[:M], phis[:M], N) err = np.sqrt(np.mean((k_series - k_recon_M)**2)) / rms_value * 100 print(f" 前{M:2d}阶谐波重建误差: {err:.4f}%")

在我这个算例里，因为刚度波形是光滑过渡的，谐波收敛非常快。前1阶谐波重建就能抓住波形的大致轮廓，前5阶时相对误差已经小于1%，前10阶时误差基本可以忽略。

但要注意，如果刚度曲线来自有限元结果，包含大量的微观波动和接触非线性，谐波收敛速度会慢很多，可能需要取到30阶甚至更高才能满足精度要求。所以程序里保留谐波阶次这个参数，让用户自己根据误差收敛情况来定，是很必要的。

3.4 典型算例结果与工程解读

我用上面的参数跑了一组实际结果，整理成下面这张表：

从表里可以读出两个重要结论：

第一，第一阶谐波幅值占绝对主导，占比超过70%。这意味着齿轮箱在啮合频率fm处的振动响应会明显高于倍频处，这和试验测得的阶次谱特征完全吻合。做减振降噪设计时，优先控制第一阶谐波的幅值，是效率最高的路径。

第二，随着阶次升高，幅值快速衰减。前5阶谐波已经囊括了超过95%的能量，前10阶之后的高阶成分基本可以忽略。这给数值仿真提供了一个重要依据：做后续动力学响应计算时，时变刚度取前5到10阶谐波就够了，不需要把整条曲线都塞进求解器，大幅提高计算效率。

4. 程序实现中的高频坑与实战经验

4.1 频谱泄漏：整周期截断是底线

某个项目里我调试一个齿轮修形方案的刚度程序，把有限元导出的刚度曲线直接丢进FFT，跑出来的谐波系数忽大忽小，完全没有规律。排查了半天，发现是有限元模型里啮合起始位置和终止位置没对齐，截取出来的序列首尾不连续，相当于在周期信号上硬切了一刀，频谱泄漏得一塌糊涂。

解决办法有两个：

第一个办法是强制整周期。输出刚度结果前，确认齿轮副旋转角度严格跨过一个基节距，也就是一个啮合周期。首尾两点理论上应该重合，但数值上由于有限元计算误差可能有微小偏差，最好做一次线性修正，强制k[N]取k[0]的值。

第二个办法是做平滑拼接。如果确实无法精确控制截断位置，可以在序列两端各取一小段做窗函数渐变处理（比如Hann窗的上升沿和下降沿），强制首尾趋近于同一个值。这样虽然会损失一点真实波形信息，但能保住频谱的干净度。实际工程中，我更推荐第一种办法，信息损失最小。

4.2 采样点数与谐波阶次上限的关系

有些朋友问我，为什么同一个模型，用512个采样点算出来的高阶谐波和用2048个点算出来的对不上。答案是混叠。

根据奈奎斯特采样定理，采样点数为N时，FFT能分辨的最高频率是N/2对应的谐波阶次。所以：

• N=512时，最多只能可靠提取到第255阶谐波；
• N=2048时，最多可以提取到第1023阶谐波。

但这只是理论上限。实际中，采样点在过渡区的分辨率直接决定了高频成分的捕获能力。如果你关心的最高谐波阶次是M，我建议采样点数至少取20M以上，也就是每个最小波形特征里至少有10个采样点，才能保证幅值精度。以我的经验，N取1024或2048对绝大多数齿轮刚度分析场景都够用，不需要无脑加大。

4.3 归一化的方向性错误

FFT归一化错误我见过太多次，这里再强调一遍。DFT的绝对值除以N，才是该频率分量在时域中的真实幅值。如果你用了2/N（常见于功率谱的处理），那确实也得到幅值，但只限于正频率单边谱的情况，此时负频率的能量已经叠加进来了。

做刚度谐波提取时，我统一的处理习惯是：

K = np.fft.fft(k_series) / N # 双边归一化 amp_dc = np.real(K[0]) # 直流分量 amp_m = 2 * np.abs(K[m]) # 第m阶谐波幅值（m>=1）

这个组合在实序列下才是物理意义上正确的幅值。你可以用正弦函数验证：如果输入是幅值为A的余弦波，按照上述方式提取出来的amp_m必然等于A，误差在浮点精度范围内。

4.4 直流分量与波动分量要区别对待

还有一个容易被忽视的问题：平均刚度k₀和波动分量，在动力学分析中的角色完全不同。

k₀决定的是齿轮副的固有频率。做模态分析或者计算临界转速时，用k₀代入即可。波动分量则是激励项，每一阶谐波幅值对应一个激励频率和幅值。如果把k₀也当成激励参与计算，会在0Hz处引入一个虚假的大幅值激励，导致系统响应计算结果完全错误。

我的程序里会单独输出一行平均刚度，并在谐波表开头注明“以下为激励分量，不含直流项”，避免下游分析时用错。

4.5 重合度接近整数时的特殊现象

重合度刚好等于整数时，刚度波形的某些谐波会异常小。比如重合度等于2.0时，理论上双齿啮合区占比100%，刚度恒定，所有谐波幅值都趋近于零。实际齿轮设计中重合度很少精确等于整数，但1.95、2.05这种接近整数的值经常碰到，这时候某些阶次的谐波幅值会非常小，形成频谱上的“凹陷”。

这个现象工程上很有价值。通过微调齿数或变位系数让重合度靠近整数，可以有效压低某个特定阶次谐波的激励，这就是齿轮修形和参数优化降噪的基本原理之一。我见过有人专门写优化程序，把重合度作为优化变量，目标函数就是第一阶谐波幅值最小化，效果非常显著，齿轮箱噪声能压下去3到5个分贝。

5. 傅立叶展开程序的实际应用延伸

5.1 动力学响应计算的输入参数化

拿到谐波系数后，一个典型的应用是齿轮系统扭转振动响应计算。把时变刚度写成：

k(t) = k₀ + Σ|K_m|·cos(2π·m·fm·t + φ_m)

代入系统的运动微分方程，用谐波平衡法或者Newmark-β直接数值积分，就能得到系统的动态响应。谐波展开的好处是，你可以很直观地看出每一阶谐波对响应的贡献，定位哪个阶次的激励和系统固有频率接近，从而引发共振。这时候针对那一阶谐波做优化，比盲目改参数高效得多。

5.2 齿轮修形参数的多目标优化

另一个重要应用是修形参数优化。齿廓修形量、修形长度等参数直接影响刚度波形，进而影响各阶谐波幅值。把傅立叶展开程序封装成一个子函数，外面套一层优化算法（遗传算法、粒子群算法都行），目标函数可以同时包含第一阶谐波幅值、整体波动幅值和传递误差，约束条件是修形量不能导致齿根弯曲强度下降过多。这种优化跑一轮，能帮你在设计阶段就避掉很多NVH风险，省下后期反复做样件测试的时间和费用。

5.3 与有限元数据的联动

我现在的标准做法是：方案阶段用解析法+傅立叶展开快速筛选参数组合，选出若干个优秀的候选方案；再用有限元细算验证这少数几个方案的刚度曲线；最后把有限元刚度数据导入傅立叶展开程序，做高精度的谐波分析和对标。这个流程兼顾了速度和精度，实际项目里跑下来非常顺手。

说到底，齿轮啮合刚度的傅立叶级数展开程序不复杂，但它是连接“几何设计”和“动力学性能”的桥梁。把这条链路跑通，齿轮箱的很多振动噪声问题，在设计阶段就能提前预判和规避。

我个人的体会是，这类程序最大的价值不在于代码本身有多巧妙，而在于你是否真正理解了每一阶谐波背后的物理含义。当你看到第1阶谐波幅值偏大时，能立刻联想到单双齿交替的过渡过程是否太突兀，能想到去调整重合度或者修形参数；当你看到某个倍频处响应超标时，能回头去查对应阶次谐波的相位是否和结构模态发生了耦合。这种“从频谱回到设计”的闭环能力，才是傅立叶级数展开程序真正值得花时间去研究的地方。

最后分享一个实用小技巧：写完傅立叶展开模块后，先用一个标准的正弦波和方波做单元测试，确认幅值和相位提取完全正确，再接入真实的刚度数据。这样可以把算法本身的问题和数据问题隔离开来，排查起来会快得多。这个习惯帮我节省了大量调试时间，建议你也在项目里用起来。