企业官网建设流程全解析

简介：面向机械工程与故障诊断领域研究人员和技术人员，这份资料系统梳理齿轮磨损故障的动态响应特征与诊断指标构建方法。内容以docx文档形式呈现，共1个文件、压缩包仅60KB，却涵盖从理论建模到代码实现的完整闭环：基于Archard磨损模型计算齿面磨损深度，采用势能法求取时变啮合刚度，将磨损等效为齿廓偏差，进而建立齿轮传动动力学模型，定量揭示磨损对啮合频率及其谐波、转频边频带的影响规律；在此基础上提出四个面向振动信号啮合频率边带的诊断指标，并验证其有效性与鲁棒性。附录包含完整Python仿真代码与逐步解释，可直接用于齿轮箱状态监测系统的设计与验证。目前已有74人学习，适合希望掌握齿轮磨损机理分析、振动信号处理及诊断指标应用的研究者参考。

1. 齿轮磨损故障动态响应特征：为什么看边带比看幅值更靠谱

齿轮磨损故障动态响应特征这条线，我在仿真正式跑通之前其实走了不少弯路。最早拿到振动数据只盯啮合频率峰值，结果发现载荷一变，幅值就飘，磨没磨损根本说不清。直到把论文里的思路落成代码才明白：磨损真正留下痕迹的地方，是啮合频率两侧的边频带。这份资源不是给你一个黑匣子，而是把从Archard磨损模型、势能法时变啮合刚度、六自由度集中参数动力学模型，到四个诊断指标计算的完整链路都拆开了，每一步都有可运行的Python代码。适合两类人：一是做齿轮箱状态监测、故障诊断算法预研的工程师，想看机理怎么变成特征；二是拿这个方向做课题或毕业设计的学生，需要一条能从理论推导到仿真验证的完整主线。它的价值不在于代码本身多复杂，而在于把“磨损→刚度变化→振动响应变化→频谱特征→诊断指标”这条因果链完整打通。

2. 磨损深度与时变啮合刚度：Archard模型和势能法的代码落地

2.1 齿轮参数类：先把工况基准算清楚

论文代码的第一步不是建模型，而是把齿轮的几何参数和工况参数集中到一个类里。这个做法很工程化，后面所有函数都从这个参数对象取值，改工况时只动一处，不会满代码找魔数。

class GearParameters: def __init__(self): # 齿轮1参数 self.z1 = 30 # 齿数 self.m1 = 2e-3 # 模数(m) self.b1 = 10e-3 # 齿宽(m) self.E1 = 2.06e11 # 弹性模量(Pa) self.v1 = 0.3 # 泊松比 # 齿轮2参数 self.z2 = 50 self.m2 = 2e-3 self.b2 = 10e-3 self.E2 = 2.06e11 self.v2 = 0.3 # 通用参数 self.pressure_angle = np.deg2rad(20) # 压力角(rad) self.H = 2.5e9 # 材料硬度(Pa) self.k = 2e-8 # 磨损系数 self.T1 = 50 # 输入扭矩(Nm) self.rpm = 1800 # 转速(rpm) # 派生参数 self.r1 = self.m1 * self.z1 / 2 # 分度圆半径(m) self.r2 = self.m2 * self.z2 / 2 self.omega1 = self.rpm * 2 * np.pi / 60 # 角速度(rad/s) self.omega2 = self.omega1 * self.z1 / self.z2 self.mesh_freq = self.rpm * self.z1 / 60 # 啮合频率(Hz)

这个类里最关键的不是齿数模数，而是三个派生参数：分度圆半径、角速度、啮合频率。它们决定了后面整个频域分析的分辨率基准。按默认参数算，输入轴转速1800rpm，齿数30，啮合频率是900Hz；齿轮1的轴转频是30Hz。这两个数字要刻在脑子里，因为后面所有边带搜索、谱线定位都是围绕它们展开的。

我一般拿到真实齿轮箱的第一件事，是把齿数和转速换成实际值重算啮合频率。这里有个容易出错的点：啮合频率的单位是Hz，不是rad/s，代码里用的是rpm * z1 / 60，如果直接用omega1去算就会差2π倍，后面边带全找不到。

2.2 Archard磨损模型：磨损深度的量级估算

Archard模型是整个磨损仿真的起点，它把磨损深度和接触压力、滑动距离、循环次数建立了定量关系。论文里用它算齿面磨损深度，而后面的刚度计算、动力学仿真、诊断指标全部以此为输入。

def archard_wear_model(params, contact_pressure, sliding_distance, cycles): """ Archard磨损模型计算磨损深度 参数: params: 齿轮参数对象 contact_pressure: 接触压力(Pa) sliding_distance: 滑动距离(m) cycles: 循环次数 返回: 磨损深度(m) """ wear_depth = params.k * contact_pressure * sliding_distance * cycles / params.H return wear_depth

公式本身很朴素：磨损深度 = 磨损系数 × 接触压力 × 滑动距离 × 循环次数 / 材料硬度。物理含义是，磨损量和摩擦功成正比，材料越硬越耐磨。默认参数下，接触压力取1e8Pa量级、滑动距离取1e-4m、循环1e5次，磨损深度大约是8微米左右，正好落在代码仿真设置的5~10微米区间内，说明参数自洽。

实际工程里有两个地方需要自己标定：一是磨损系数k，齿轮材料、润滑状态、表面粗糙度都会影响它，论文代码给的2e-8是典型的干摩擦经验值，如果是油润滑，通常要往下调一个量级；二是滑动距离，齿轮啮合过程中齿面各点的滑动速度不一样，节线附近滑动距离趋近于零，齿根齿顶最大，所以真正的磨损分布不均匀，但做诊断指标研究时先取等效值就够了。论文后面的扩展代码里把接触点沿齿廓分布成100个点，计算磨损分布而不是单点磨损深度，这个扩展在磨损形态分析时有用，做状态监测综合指标时不必上。

2.3 势能法时变啮合刚度：简化版与全量版

时变啮合刚度是齿轮动力学里最核心的时变参数，它随啮合位置周期性变化，正是这个周期性变化在频域里制造出啮合频率及其谐波。论文里用势能法计算，原理解释一下：啮合轮齿可以看成悬臂梁，承受接触载荷时产生弯曲变形、剪切变形、轴向压缩变形和局部接触变形，每一部分都对应一个刚度分量，最后串联成总啮合刚度。

def time_varying_mesh_stiffness(params, theta, wear_depth=0): """ 计算时变啮合刚度 参数: params: 齿轮参数对象 theta: 齿轮1转角(rad) wear_depth: 磨损深度(m) 返回: 啮合刚度(N/m) """ contact_ratio = 1.5 # 假设重合度 mesh_period = 2 * np.pi / params.z1 # 啮合周期 base_stiffness = 1e8 # 基础刚度(N/m) stiffness_variation = 0.2 * np.sin(params.z1 * theta) # 刚度波动 wear_effect = 0.5 * (wear_depth / 1e-5) if wear_depth > 0 else 0 total_stiffness = base_stiffness * (1 + stiffness_variation - wear_effect) return total_stiffness

这个简化版把时变啮合刚度拆成三项：常数项1e8、波动项0.2倍的啮频正弦、磨损项。其中波动项模拟的是啮合过程中同时啮合的齿对数变化——单双齿交替导致刚度周期性起伏，这是啮频激励的主要来源。磨损项让刚度随磨损深度线性下降，每10微米磨损对应刚度下降50%，这个系数是为了让仿真在5~10微米磨损下产生可观察的谱变化。

提示：磨损对刚度的实际影响不是纯线性的，但状态监测场景下，早期磨损阶段用线性近似完全够用。真要计算每个齿的精确刚度，得用有限元或解析公式逐齿廓积分，工程上一般到后期故障分析才做。

论文扩展部分给的MeshStiffnessCalculator是更完整的势能法实现，把赫兹接触刚度、弯曲刚度、轴向压缩刚度、剪切刚度分别算出来再串联：

# 完整势能法的四项刚度（代码扩展部分） E_eq = 2 / ((1 - v1**2) / E1 + (1 - v2**2) / E2) # 等效弹性模量 k_h = np.pi * E_eq * b1 / (4 * (1 - v1**2)) # 赫兹接触刚度 L = r1 * theta # 接触点到齿根距离 k_b = b1 * E1 / (6 * L**3) # 弯曲刚度 k_a = b1 * E1 / L # 轴向压缩刚度 k_s = b1 * E1 / (2 * (1 + v1) * L) # 剪切刚度 k_total = 1 / (1/k_h + 1/k_b + 1/k_a + 1/k_s) # 串联

四项刚度是串联不是并联，因为轮齿变形是各分量变形的叠加，柔度相加再取倒数。弯曲刚度占主导，因为它对接触点到齿根距离L的三次方敏感。做工程诊断时如果只需要啮频边带特征，简化版的正弦波动就够；如果要研究磨损分布对谐波成分的精细影响，建议用完整版。

3. 动态响应仿真：集中参数模型与odeint的六自由度实现

3.1 状态变量与运动方程：为什么是12维

齿轮系统的动力学模型采用集中参数法，把每个齿轮简化为质量和转动惯量，同时考虑水平和垂直方向的平移自由度，以及绕轴转动自由度。论文代码里把两个齿轮完全独立建模，每个齿轮3个自由度，共6个自由度，对应12个状态变量：6个位移（x、y、theta）和6个速度。

def gear_system_equations(y, t, params, wear_depth): # 解包状态变量 x1, y1, th1, x2, y2, th2, dx1, dy1, dth1, dx2, dy2, dth2 = y # 当前啮合刚度和静态传递误差 km = time_varying_mesh_stiffness(params, th1, wear_depth) STE = static_transmission_error(params, th1, wear_depth) # 啮合线上相对位移（含传递误差激励） delta_th = params.r2 * th2 - params.r1 * th1 - STE # 啮合力 Fm = km * delta_th # 动力学方程 ddx1 = Fm * np.cos(params.pressure_angle) / m1 ddy1 = Fm * np.sin(params.pressure_angle) / m1 - 9.8 ddth1 = (T1 - Fm * params.r1) / I1 ddx2 = -Fm * np.cos(params.pressure_angle) / m2 ddy2 = -Fm * np.sin(params.pressure_angle) / m2 - 9.8 ddth2 = (T2 + Fm * params.r2) / I2

最核心的一行是delta_th = params.r2 * th2 - params.r1 * th1 - STE，它表示两齿轮在啮合线方向上的相对位移，减去静态传递误差后才是弹性变形量。啮合力等于时变刚度乘以这个弹性变形量，这就是磨损影响振动的传导路径：磨损改变了STE和刚度，继而改变啮合力，最后改变振动响应。这里需要理解：STE是位移激励，刚度变化是参数激励，两者同时作用在系统上。

初始条件里两个角速度直接设为稳态值，这个细节很关键。如果从零转速启动，系统会经历剧烈的瞬态冲击，频谱里全是启动暂态成分，掩盖磨损特征。工程上做齿轮箱振动分析也一样，要等转速稳定后再取数据。

3.2 odeint求解：步长、时长与频率分辨率

动态响应求解用scipy的odeint，这是求解常微分方程最常用的函数，内部自动变步长，只需要给定初始条件、时间数组和传入参数。

t = np.arange(0, t_end, dt) initial_conditions = [0, 0, 0, 0, 0, 0, 0, 0, params.omega1, 0, 0, params.omega2] sol = odeint(gear_system_equations, initial_conditions, t, args=(params, wear_depth)) theta1 = sol[:, 2] theta2 = sol[:, 5] TE = params.r2 * theta2 - params.r1 * theta1

时间数组设置是仿真成败的关键参数。t_end=0.1和dt=1e-5的组合带来两个结果：一是求解步数10000步，odeint处理这个规模很轻松；二是FFT的频率分辨率等于1/t_end，也就是10Hz。默认转频30Hz，啮频900Hz，10Hz分辨率意味着每条边带之间隔3根谱线，勉强能看到边带结构。

实际复现时我建议把t_end加到0.5秒，分辨率降到2Hz，边带结构会清晰很多。代价是求解时间变长，但对单次仿真来说完全可接受。dt取1e-5是数值稳定性的考量，如果缩小到1e-6会更稳但计算量翻倍，如果放大到1e-4，odeint可能报错。参数调整原则就一句话：先保数值稳定，再提频率分辨率。

3.3 时域与频域对照：磨损在频谱上留下什么

仿真输出后，论文代码用FFT分别绘制不同磨损深度下的时域和频域响应。时域看的是传递误差波形，频域看的是啮频成分与边带结构的变化。磨损的影响不是让整体幅值均匀变化，而是集中在啮合频率及其谐波两侧。

幅值调制机理值得展开：磨损导致齿面形貌改变，啮合刚度每一转出现一次局部变化。啮合频率是载波频率900Hz，轴转频是调制频率30Hz，调制结果就是在900Hz两侧对称出现30Hz间隔的边频带。磨损越重，这个局部变化的强度越大，边带能量越突出。论文的结论“磨损主要影响啮合频率及其谐波成分，边频带以转频为主且随磨损增加而变化”，就是从这张频谱图里提炼出来的。

注意：如果仿真参数里齿轮2是输出轮，它的转频是18Hz，边带间隔会随之改变。判断边带属于哪个齿轮，先算两个轴的转频，再对照谱线间隔，不要默认边带一定等于输入轴转频。

4. 四个诊断指标：边带特征如何变成可计算的数

4.1 指标定义与选择依据

论文的核心贡献之一是把频谱边带特征压缩成四个可用数值比较的诊断指标。为什么要压缩？因为直接对比频谱曲线在工程上不可行——振动幅值受载荷、转速、传感器灵敏度影响，不同工况下的频谱不能直接相减。指标的设计原则是：幅值类指标反映磨损绝对强度，比值类指标消除工况影响。

选这四个指标的逻辑是：I1和I3盯住载波能量，I2和I4盯住调制深度。工程上我更看重I2和I4，因为比值天然消除了传感器增益和载荷波动的影响。论文验证了这组指标在不同磨损状态下的有效性和鲁棒性，实际做状态监测阈值报警时，边带比类指标比绝对幅值稳定得多。

4.2 实现细节：边带范围和频率分辨率

指标计算的代码不长，但有几个细节决定结果靠不靠谱。

def calculate_wear_indicators(spectrum, freq, mesh_freq, shaft_freq): mesh_idx = np.argmin(np.abs(freq - mesh_freq)) mesh2_idx = np.argmin(np.abs(freq - 2 * mesh_freq)) I1 = spectrum[mesh_idx] # 边带搜索范围：以转频为半宽 sideband_range = int(shaft_freq / (freq[1] - freq[0])) left_side = np.max(spectrum[mesh_idx - sideband_range:mesh_idx]) right_side = np.max(spectrum[mesh_idx:mesh_idx + sideband_range]) I2 = (left_side + right_side) / I1 return I1, I2, I3, I4

sideband_range的计算逻辑是用转频除以频率分辨率，得到边带搜索窗口覆盖多少根谱线。默认参数下shaft_freq=30Hz，分辨率=10Hz，窗口就是3根谱线。这里有个隐患：如果FFT频率分辨率太粗，转频可能落在谱线间隙里，实际边带峰值被相邻谱线均摊，np.max取到的值偏小。我的经验是边带搜索前先用插值或加窗处理频谱，至少要让窗口内包含5根以上谱线，否则I2的数值波动会很大。

4.3 主程序里那个“演示频谱”是怎么回事

原代码主程序里计算指标用的spectrum不是仿真输出的FFT结果，而是手工构造的高斯函数近似频谱。这是为了单独演示指标函数的功能，但很容易让人误以为仿真和指标已经自动串联。复现时的关键一步是把两段接起来：把simulate_gear_dynamics里计算出的FFT幅值保存为变量，传入calculate_wear_indicators。

# 仿真后取频谱幅值 yf = fft(TE) xf = fftfreq(n, dt)[:n//2] spectrum = 2.0 / n * np.abs(yf[:n//2]) # 注意去除直流分量 spectrum = spectrum[xf > 0] xf = xf[xf > 0] # 再计算四个指标 I1, I2, I3, I4 = calculate_wear_indicators(spectrum, xf, params.mesh_freq, 30.0)

接上之后，整个链路才真正闭环：磨损深度 → 刚度/STE变化 → 动力学响应 → FFT频谱 → 诊断指标。这个闭环是做磨损趋势分析的基础，后面验证指标单调性时也必须走这条路。

5. 复现排查：频谱读数异常与仿真发散的处理记录

5.1 零频处出现巨大峰值，量程被压扁

现象：FFT频谱画出来后，0Hz附近有一个极大的尖峰，其他频段全部贴地，啮频和边带完全看不清。

原因：传递误差TE存在直流分量。初始条件、重力项和系统静变形都会在时域响应中叠加一个常数偏置，FFT对直流分量极其敏感，幅值远高于交流成分。

解决：FFT之前先对TE序列去均值。一行代码就行：

TE = TE - np.mean(TE) # 去直流后再做FFT

去均值之后，频谱纵轴自动缩放，啮频和边带结构就出来了。做现场振动分析时同理，加速度信号进入FFT前先做高通滤波或去趋势，处理的是同一个问题。

5.2 odeint求解发散，输出变成inf或NaN

现象：仿真运行过程没有报错，但频谱图里数值异常，检查TE数组发现中间一段直接变成无穷大。

原因：这是数值刚性问题，啮合刚度1e8量级，而阻尼系数cm=1e3相对太小，系统等效为高刚度低阻尼，微分方程的特征值跨度大，odeint默认的求解器在局部区间步长控制不住。另外初始角速度如果和稳态不匹配，系统一开始就有很大的瞬态冲击。

解决：两条路都试过，最有效的是把阻尼系数从1e3提到1e4，增加能量耗散；如果还发散，换scipy.integrate.solve_ivp并指定method='Radau'，Radau方法对刚性方程比odeint默认的LSODA更稳。另外检查初始条件里omega1和omega2是否满足传动比关系，不匹配等于给了系统一个初始冲击。

5.3 频谱边带间隔对不上转频

现象：频率图上900Hz两侧确实有边带，但间隔是60Hz而不是预期的30Hz。

原因：转频算错。边带间隔等于轴转频，角速度omega1除以2π才是转频。如果直接把omega1的数值188.5当转频用，或者把转速1800rpm除以60忘记换算，间隔都会差一倍。另一个隐蔽原因是边带搜索窗口设成了2倍转频宽，取max时把相邻两个边带都包进去了。

解决：用f = omega1 / (2 * np.pi) = 30Hz作为基准，先打印确认再计算sideband_range。我现在的习惯是在参数类里直接加一个shaft_freq属性，避免每次手动换算。

5.4 不同磨损深度下频谱几乎没变化

现象：把wear_depths设为[0, 5e-6, 10e-6]跑仿真，三条频谱曲线几乎重合。

原因：磨损项对刚度的修正量没起作用。检查时发现time_varying_mesh_stiffness里wear_effect的计算依赖wear_depth传入，如果主程序调用时忘了传参，wear_depth一直为0，磨损效应完全不生效。另一个可能是磨损深度量级太小，5微米对1e8的刚度来说只改变0.25%，频谱差异确实不明显。

解决：先确认wear_depth确实传进了刚度函数，打印一下刚度数值验证磨损项在起作用。然后把磨损深度序列拉大，比如[0, 5e-6, 20e-6, 50e-6]，先验证趋势，再回到实际量级。

5.5 指标计算和仿真频谱各算各的

现象：四个指标打印出来数值一直不变，改磨损深度也影响不了它。

原因：主程序最后计算指标用的spectrum是手工构造的示例频谱，跟仿真结果完全没关系。指标模块没问题，但它没有接在动力学仿真后面，所以磨损深度怎么改都影响不到指标数值。

解决：按4.3节的方式把仿真FFT结果保存下来，传入calculate_wear_indicators。这条是我复现时最有感触的地方——论文代码为了分模块演示，把链路的最后一环切断了，如果不主动接上，很容易在演示阶段就以为任务完成了。

6. 验证诊断指标：单调性检验与台架试验对照

6.1 单调性检验：指标要随磨损深度稳步上升

验证指标有没有用，第一件事是做单调性检验。把磨损深度序列打密，例如[0, 2e-6, 5e-6, 10e-6, 20e-6, 40e-6]，对每个深度跑一遍仿真并计算I2和I4，观察是否随磨损单调变化。

wear_levels = [0, 2e-6, 5e-6, 10e-6, 20e-6, 40e-6] indicator_results = [] for wd in wear_levels: sol = odeint(gear_system_equations, initial_conditions, t, args=(params, wd)) # 提取TE、去均值、FFT TE = params.r2 * sol[:, 5] - params.r1 * sol[:, 2] TE = TE - np.mean(TE) spectrum = 2.0 / n * np.abs(fft(TE)[:n//2]) I1, I2, I3, I4 = calculate_wear_indicators( spectrum, xf, params.mesh_freq, 30.0) indicator_results.append((I2, I4))

判断标准我一般定两条：一是I2、I4随磨损深度严格单调上升，允许中间有小波动但不允许下降后无法恢复；二是重复跑三次仿真，指标值的波动不超过15%。满足这两条，指标才敢用到监测系统里做趋势报警。

6.2 与台架试验对照：频谱特征位置必须一一对上

仿真通过了还要过试验这一关。加速度传感器装在齿轮箱轴承座垂直方向，采样率至少8kHz，覆盖到4倍啮合频率以上，防止谐波混叠。试验数据先算转速，精确到0.1Hz，再算啮合频率；FFT参数和仿真保持一致，分辨率不低于10Hz。对照时主要看三处：啮合频率处有无明显峰值、边带间隔是否等于轴转频、磨损加剧后边带幅值是否相对啮频上升。有一个对不上就回到仿真侧查参数，通常问题出在转速波动——试验时转速不稳会把边带抹平，这是现场和仿真最大的差别。

6.3 我的验证习惯

这套流程走下来，我形成了一个固定习惯：每次跑齿轮动力学仿真，第一件事就是算频率分辨率Δf和转频的比值，确认Δf小于转频的一半，否则边带峰值会被相邻谱线稀释，指标计算必然失真。这个习惯帮我避掉了很多次“数据没问题但指标不好看”的情况。从那以后我每次做频谱相关的诊断分析，都强制走一遍这个检查。希望帮到你。

本文还有配套的精品资源，点击获取