搞齿轮故障诊断这些年经常被同行问到的一句话是剥落故障在振动信号上到底长什么样说实话光靠现场测振谱很难讲清楚因为传感器采到的响应是齿轮副时变啮合刚度、齿侧间隙、误差激励、负载波动等多因素耦合后的综合结果。我一直建议做故障诊断的工程师把“含剥落故障直齿轮啮合刚度”的计算先写成程序再把刚度结果喂进齿轮非线性动力学方程里做数值仿真这样才能把故障特征从一堆耦合信号里单独拧出来看清楚。本文就从啮合刚度的机理推导开始讲到非线性动力学模型的搭建和MATLAB程序实现最后附上我实际调试中踩过的坑希望能给正在做齿轮动力学仿真的朋友一些参考。这篇文章适合三类人一是做齿轮箱故障诊断、想在仿真里复现剥落特征的研究生和工程师二是做传动系统设计、需要评估损伤对动态性能影响的机械工程师三是刚接触齿轮非线性动力学、想把论文里的方程变成能跑出结果的程序的初学者。我会尽量把公式背后的物理含义讲透代码结构也会拆开说明。1. 剥落故障对啮合刚度的作用机理1.1 剥落不是“表面刮花”它直接改写齿轮副的承载路径齿轮剥落spalling属于接触疲劳失效本质是齿面次表层在循环接触应力作用下萌生裂纹裂纹扩展到一定深度后引起材料成片脱落形成不规则凹坑。很多人对剥落的理解停留在“齿面多了个坑”但从动力学仿真的角度看这个坑的存在意味着两个重要变化第一齿轮副在啮合到剥落区域时原本由整个齿宽承担的接触载荷被迫转移实际有效接触面积减小第二轮齿悬臂梁模型的截面几何发生变化弯曲刚度、剪切刚度随之改变。这两点最终都归结到一个核心物理量——啮合刚度的时变曲线出现了局部的凹坑和突变。理解这一点是后续一切工作的前提。齿轮啮合刚度并不是恒定值它随啮合位置周期变化这本身就是齿轮系统最主要的内部激励源。当剥落故障存在时相当于在原本周期变化的刚度曲线上叠加了一个“局部缺陷脉冲”这个脉冲会激励起齿轮系统的共振响应在频谱上表现为啮合频率及其谐波附近的边带。我见过不少初学者直接拿恒定刚度建模仿真出来的频谱干干净净根本复现不出边带特征就是因为没有把剥落引起的刚度变化放进去。1.2 恒定刚度假设在故障仿真里的局限经典齿轮设计中常用平均啮合刚度或ISO标准推荐的刚度公式来算承载能力这在强度校核阶段够用但做故障诊断仿真就远远不够了。齿轮副内部激励的主要来源就是时变啮合刚度TVMS, Time-Varying Mesh Stiffness恒定刚度等于把最大的激励源直接抹掉了。剥落故障诊断的核心逻辑是“故障改变了刚度→刚度改变激励→激励改变响应”如果第一步就用了恒定刚度后面全是空中楼阁。1.3 程序化解决的核心思路把剥落故障纳入动力学仿真的完整链路是建立剥落的几何参数模型 → 基于势能法计算含故障的时变啮合刚度 → 把刚度序列作为时变系数代入非线性运动微分方程 → 数值积分得到振动响应 → 对响应做频谱、相图、庞加莱截面分析。这套链路每一步都有数学依据每一步都可以用程序固化下来。我下面按这个顺序逐层展开。2. 剥落故障的数学化描述与几何建模2.1 矩形剥落模型的三要素工程中为了解析计算方便通常把不规则剥落坑简化为矩形凹槽。这个简化在故障诊断学术界是通用做法McGraw-Hill的教材和大量论文里都用这种近似。矩形剥落需要定义三个参数剥落轴向长度 (L)、剥落深度 (h_s)、剥落沿齿廓方向的位置通常用距齿顶的距离或对应滚动角描述。我以自己常用的参数为例模数 (m3)小齿轮齿数 (z_120)大齿轮齿数 (z_230)齿宽 (b20) mm压力角 (\alpha20^\circ)标准安装距。剥落参数取轴向长度 (L3) mm、深度 (h_s0.5) mm位置在小齿轮单双齿啮合交替区域附近。剥落深度不能超过齿根危险截面的合理范围否则刚度会算成负值这在第7章会专门讲。2.2 啮合过程的接触点移动轨迹要计算时变啮合刚度先得弄清楚一个啮合周期内接触点的位置怎么变化。直齿轮副啮合时接触线沿齿廓方向从齿根向齿顶或反过来扫过同时啮合的齿对数在1和2之间切换。啮合刚度的周期就是端面基节 (p_{bt} \pi m \cos\alpha) 对应的小齿轮转角 (\varphi_z 2\pi / z_1)。程序实现时我把小齿轮转过一个齿距角均分成 (N200) 个离散位置每个位置分别计算当前参与啮合的齿对数、每个齿上接触点的坐标再根据接触点是否落入剥落区间来决定是否在刚度公式中计入剥落影响几何参数。这里有个关键细节剥落位置若发生在单齿啮合区刚度下降会非常剧烈因为那时只有一对齿承担全部载荷若在双齿啮合区载荷由两对齿分担刚度下降相对温和。这个位置敏感性本身就是故障诊断里区分剥落严重程度的重要依据。2.3 剥落影响怎么进入截面参数势能法把轮齿简化成变截面悬臂梁弯曲刚度和剪切刚度都依赖于截面惯性矩和截面面积。无故障时齿廓渐开线决定了截面厚度沿高度方向的变化有剥落时在剥落位置对应的截面沿齿宽方向的有效厚度减小了 (h_s)实际等效为截面惯性矩降低。处理方式是在计算截面参数时先判断当前截面高度是否落入剥落区间再更新该截面的厚度和面积。这个判断逻辑是程序里最容易写错的地方——很多人直接把厚度减去 (h_s)却没考虑剥落只在部分齿宽上存在导致刚度下降幅度过大。3. 啮合刚度计算的势能法完整推导3.1 五种能量分量与对应的刚度表达式势能法能量法的核心思想是把一对啮合轮齿的弹性变形分解为五个部分分别计算各自对应的刚度再按照串联关系合成综合啮合刚度。五种刚度为赫兹接触刚度 (k_h)、轮齿弯曲刚度 (k_b)、剪切刚度 (k_s)、轴向压缩刚度 (k_a)、齿基体弹性刚度 (k_f)。以弯曲刚度和剪切刚度为例公式推导基于悬臂梁应变能[ U_b \int_{0}^{d} \frac{(F\cos\alpha_1(d-x) - F\sin\alpha_1 h_x)^2}{2EI_x} dx ][ U_s \int_{0}^{d} \frac{1.2 (F\cos\alpha_1)^2}{2GA_x} dx ]其中 (d) 是载荷作用点到齿根的距离(x) 是沿齿高方向的积分坐标(h_x) 是截面到中性轴的距离函数(I_x) 和 (A_x) 分别是截面惯性矩和截面积。计算单个轮齿的刚度 (k_i F^2 / (2U_i))再把五个分量的柔度刚度的倒数按串联关系相加[ \frac{1}{k_m} \frac{1}{k_h} \frac{1}{k_{b1}} \frac{1}{k_{s1}} \frac{1}{k_{a1}} \frac{1}{k_{f1}} \frac{1}{k_{b2}} \frac{1}{k_{s2}} \frac{1}{k_{a2}} \frac{1}{k_{f2}} ]下标1、2代表主从动轮。Hertz接触刚度采用 (k_h \frac{\pi E b}{4(1-\nu^2)})其中 (E) 为弹性模量(\nu) 为泊松比(b) 为齿宽。齿基体刚度 (k_f) 可以用Sainsot等的经验公式计算它反映齿圈弹性变形的影响占比不大但在计及薄腹板齿轮时不能忽略。3.2 剥落如何修正积分上下限与截面参数有剥落时积分过程中的 (I_x) 和 (A_x) 需要做局部修正。假设剥落矩形凹槽对应齿廓高度区间为 ([x_{s1}, x_{s2}])齿宽方向的剥落长度从 (0) 到 (L)则落在该区间的截面其有效齿宽从 (b) 变为 (b - L)截面厚度相应减薄。更严格的做法是用二维接触模型把剥落坑附近接触压力重新分布但解析法中常用的简化是直接修正截面参数。这虽然粗糙但在工程精度内足以捕捉到刚度曲线的下凹特征。单齿对啮合刚度计算完成后还需要根据啮合周期内单双齿啮合区的切换来合成一个完整周期的时变刚度曲线。双齿啮合区相当于两个齿对的刚度并联合成规则是[ k_{mesh}(t) \sum_{\text{同时啮合的齿对}} k_{mj}(t) ]程序里用一个双循环实现外层遍历小齿轮转角内层判断该转角下同时参与啮合的齿对编号每个齿对的接触点位置不同各自的剥落影响也不同。3.3 刚度曲线的典型特征我算过大量含剥落直齿轮副的刚度曲线典型特征非常明显在剥落对应的小齿轮转角处综合刚度值出现一个局部的“V”形下凹深度与剥落尺寸正相关宽度与剥落沿齿廓方向的长度正相关。更关键的是剥落引起的刚度损失不仅影响当前齿对还会影响相邻啮合周期的边界导致刚度曲线不再是严格周期函数这就是程序中必须把至少两个完整啮合周期的刚度算出来再截取稳态段的原因。4. 齿轮非线性动力学模型与程序实现4.1 单自由度扭转模型的建立有了刚度序列就可以建动力学方程。我采用经典的直齿轮副单自由度集中质量模型广义坐标取动态传动误差 (x r_{b1}\theta_1 - r_{b2}\theta_2 - e(t))其中 (r_b) 是基圆半径(e(t)) 是综合啮合误差。运动微分方程为[ m_e \ddot{x} c \dot{x} k(t) g(x) F_m F_a \sin(\omega t \varphi) ]其中 (m_e) 是等效质量(c 2\zeta \sqrt{m_e \bar{k}}) 是啮合阻尼(\zeta) 取0.01~0.05(\bar{k}) 是平均啮合刚度(F_m) 是平均载荷(F_a) 是载荷波动幅值(\omega) 是激励频率。(g(x)) 是齿侧间隙函数[ g(x) \begin{cases} x - b_g, x b_g \ 0, -b_g \le x \le b_g \ x b_g, x -b_g \end{cases} ]其中 (b_g) 是半齿侧间隙。这个非线性项直接导致齿轮系统出现跳变、混沌等复杂动力学行为是“非线性动力学”的根源所在。方程两边的量级差别很大直接积分会碰到数值困难。我通常先做无量纲化令 (x_n x / b_c)(b_c) 为特征长度取间隙量级(\tau \omega_n t)(\omega_n \sqrt{\bar{k}/m_e})。无量纲化之后方程变成[ x_n 2\zeta x_n \frac{k(\tau)}{\bar{k}} g(x_n) \frac{F_m}{m_e b_c \omega_n^2} \frac{F_a}{m_e b_c \omega_n^2} \sin(\Omega \tau) ]这里 (g(x_n)) 也要同步除以 (b_c)。无量纲化有两个好处一是把数值量级统一到 (10^{-1}) 到 (10^1) 之间显著降低积分器的绝对误差控制难度二是让结果具有通用性便于以后换参数时做无量纲对比。这一步我强烈建议不要偷懒跳过直接拿SI单位积分时位移量级在 (10^{-6}) m 左右刚度量级在 (10^8) N/m两者跨了14个数量级ode45很容易把时间步长压到极小导致计算时间爆炸。4.2 时变刚度与误差激励的程序化处理程序里时变啮合刚度不能写成解析表达式因为剥落导致它不光滑。我采用查表法先在预处理阶段计算出两个完整啮合周期的刚度序列 (k_m[i])(i1,...,N)(N400)通过线性插值得到任意时刻的刚度值。误差激励 (e(t)) 包含齿频误差和转频误差通常用简谐波叠加模拟[ e(t) e_0 e_{hf}\sin(2\pi f_m t \phi_1) e_{lf}\sin(2\pi f_r t \phi_2) ]其中 (f_m) 是啮合频率(f_r) 是转频(e_0) 是常值误差(e_{hf}) 和 (e_{lf}) 分别是高低频误差幅值。程序里把这些都放进全局参数结构体方便批量修改。4.3 积分器的选择与参数设置微分方程数值积分我用MATLAB的ode45但对含间隙的非线性系统ode45在某些刚度过大突变的位置会触发变步长算法的“事件检测”机制导致步长剧烈震荡。如果碰到这种情况有三个办法一是改用ode23s针对刚性方程二是把刚度序列做轻度的平滑滤波消除数值计算的微小跳变三是明确规定RelTol和AbsTol我一般设RelTol1e-6, AbsTol1e-7。实测下来对单自由度齿轮系统ode45配合适度平滑的刚度输入就足够了ode23s反而在稳态段耗时更长。积分时长方面至少需要让系统跑完200个啮合周期再丢弃前50个周期的瞬态响应取后150个周期做分析。很多人只跑几十个周期就拿来分析频谱结果瞬态分量混在频谱里边带特征一团模糊。这个时长问题在程序里直接用一个循环控制总积分时间t_end 250 * T_mesh后处理时从50 * T_mesh开始截取数据。5. 程序架构与核心代码实现5.1 模块划分与数据流完整的程序我拆成四个文件职责清晰主程序main_gear_dyn.m负责定义参数、调用各模块、绘图刚度计算函数compute_mesh_stiffness.m只做一件事——输入小齿轮转角序列和剥落参数输出对应时刻的啮合刚度序列动力学右端项函数gear_ode_right.m接收状态变量和时间返回导数后处理脚本plot_results.m专门画频谱、相图、庞加莱截面。模块划分的好处是换一组参数、换一种剥落尺寸、换一个间隙值时不需要改动核心逻辑只改参数区就行。5.2 刚度计算函数核心逻辑先看刚度计算函数的框架function k_mesh compute_mesh_stiffness(gear_params, spall_params, theta_seq) % 输入 % gear_params: 齿轮几何与材料参数结构体 % spall_params: 剥落参数结构体长度L、深度h_s、角位置theta_sp % theta_seq: 小齿轮转角序列1×N % 输出 % k_mesh: 综合啮合刚度序列1×N N length(theta_seq); z1 gear_params.z1; z2 gear_params.z2; theta_pitch 2*pi/z1; % 啮合周期对应的转角 k_mesh zeros(1, N); for i 1:N theta theta_seq(i); % 计算当前转角对应的基节内位置 phi mod(theta, theta_pitch) / theta_pitch; % 0~1 归一化啮合位置 % 判断是单齿啮合还是双齿啮合区 [pair_count, contact_params] detect_contact_zone(phi, gear_params); k_pair_sum 0; for j 1:pair_count % 对每个接触齿对计算含剥落修正的单齿对刚度 k_pair_sum k_pair_sum single_pair_stiffness(contact_params(j), ... gear_params, spall_params); end k_mesh(i) k_pair_sum; end end这个函数里最容易出错的是detect_contact_zone的判断逻辑。直齿轮啮合重合度通常在1到2之间也就是说大部分时间有两对齿同时啮合只有一小段是单齿啮合。判断依据是啮合位置相对于基节 (p_{bt}) 的比值基节范围内前半段是双齿啮合中间是单齿啮合后半段又是双齿。剥落修正的核心在single_pair_stiffness里。它内部会先计算当前接触点对应的齿廓高度坐标判断该高度是否落在剥落区间 ([h_{s1}, h_{s2}])然后决定积分时的截面参数function k_single single_pair_stiffness(cp, gear_params, spall_params) % cp: 接触点几何信息含d, x范围等 % 无剥落时的截面参数 I_x calc_inertia(cp, gear_params); % 截面惯性矩 A_x calc_area(cp, gear_params); % 截面面积 if is_in_spall_zone(cp, spall_params) % 关键判断 b_eff gear_params.b - spall_params.L; % 有效齿宽 I_x I_x * (b_eff / gear_params.b); A_x A_x * (b_eff / gear_params.b); % 深度修正等效厚度减薄 I_x I_x * (1 - spall_params.h_s / cp.h_sec)^3; % 按矩形截面厚度立方修正 end % 计算各项应变能并累加... end这里有个细节矩形截面惯性矩与厚度立方成正比所以深度修正用了三次方比例这是很多初学程序最容易漏掉的地方漏掉后剥落对刚度的影响会被严重低估。5.3 动力学方程右端项与主循环右端项函数的核心代码如下function dydt gear_ode_right(t, y, p) % y(1)x_n, y(2)x_n % p: 结构体含k_seq, theta_seq, gap等 theta mod(p.omega_n * t, 2*pi/p.z1); k_now interp1(p.theta_seq, p.k_seq, theta, linear); % 间隙函数 if y(1) p.bg_n g_val y(1) - p.bg_n; elseif y(1) -p.bg_n g_val y(1) p.bg_n; else g_val 0; end dydt zeros(2,1); dydt(1) y(2); dydt(2) -2*p.zeta_n*y(2) - k_now * g_val p.f_n p.fa_n*sin(p.Omega*t); end主程序只需要调用ode45并做后处理tspan [0, p.t_end]; y0 [0.1; 0]; % 初始位移和速度按无量纲量级取 [t, y] ode45((t,y) gear_ode_right(t,y,p), tspan, y0, opts); % 截取稳态段 idx find(t p.steady_start); y_steady y(idx, 1); tt_steady t(idx);6. 仿真结果分析与故障特征解读6.1 含剥落时变刚度曲线的V形下凹用我前面给的参数跑出来的刚度曲线在剥落角位置会看到明显的局部下凹凹坑深度约为正常刚度的8%~15%具体取决于剥落尺寸与齿宽的比值。剥落长度 (L) 增大时凹坑变宽变深剥落深度 (h_s) 增大时凹坑深度增长更快。这个下凹就是后续振动响应一切异常的总源头。我习惯把故障工况和健康工况的刚度曲线叠在一张图里画用阴影标注剥落区间这样演示故障机理最直观。注意刚度曲线的两端会有啮合周期切换引起的突变那是双齿啮合区与单齿啮合区交替的正常现象不要和剥落造成的下凹混淆。6.2 振动响应的时域与频域特征把剥落刚度代入动力学方程后时域位移响应的最大变化出现在剥落对应的转角附近振动幅值增大且波形出现调制包络。频域上啮合频率 (f_m z_1 n/60)(n) 为输入转速处的幅值上升同时在 (f_m \pm k f_r) 处出现明显的边带簇(k1,2,...)边带间隔等于小齿轮转频 (f_r)。边带的幅值不对称性可以用于判断剥落发生在主动轮还是从动轮上这是工程诊断里一个非常实用的判据。频谱分析程序我用标准FFT配合汉宁窗采样点数取 (2^{14})频率分辨率控制在转频的1/4以下否则边带会被谱线间隔掩盖。6.3 相图与庞加莱截面判断系统状态齿轮含间隙系统的典型非线性行为包括周期运动、拟周期运动、混沌。剥落故障改变了局部刚度可能把原本稳定的周期运动推向混沌或产生周期分岔。判断方法看相图位移-速度平面轨迹和庞加莱截面周期1运动对应相图是一条闭合曲线庞加莱截面一个映射点周期2运动相图有两条交织曲线截面两个点混沌时相图杂乱无章且永不重复截面上出现分形结构的点集。我在程序中用“每啮合周期采样一次庞加莱点”的办法记录每个 (t k T_m) 时刻的位移和速度画成散点图。剥落工况和健康工况对比时即使都处于周期运动剥落后的庞加莱映射点也会出现漂移或分裂这是早期故障检测的敏感指标。7. 常见问题与排查技巧实录7.1 刚度曲线出现负值或剧烈抖动我最初调试程序时碰到过刚度突然变成负值的情况查了两天才找到原因剥落深度 (h_s) 超过了该截面本身的厚度导致修正后的截面参数为负。齿轮齿顶和齿根过渡区的截面厚度本来就不一样越靠近齿顶截面越薄如果剥落参数设置得过大就很容易触发负刚度。解决办法是在几何参数初始化时增加一个检查函数确保所有截面在计入剥落后的有效厚度都大于判别值例如 (h_{eff} \ge 0.05 \times h_{max})。另一个导致抖动的原因是修正 (I_x) 时用了简单的if/else边界判断边界处刚度出现不连续跳变。我给剥落边界加了3~5个网格点的过渡带让截面参数平滑过渡曲线就自然了。7.2 积分不收敛或耗时过长如果ode45在剥落刚度突变处反复缩短步长先检查刚度序列是否含NaN或Inf再检查无量纲化是否正确。我调试时发现一个坑无量纲化后激励频率 (\Omega) 应该是相对啮合频率的比值很多人直接用了物理频率导致积分器在过高频率下步长被压死计算时间成倍增加。正确做法是 (\Omega \omega / \omega_n)其中 (\omega) 是啮合频率(\omega_n) 是系统固有频率。7.3 模型验证的三个层次仿真程序写完后一定要验证否则结果没有说服力。我总结了三层验证第一层在无剥落、无间隙(b_g0)的条件下刚度曲线和固有频率应与解析解或有限元结果误差在5%以内第二层把仿真位移响应的平均值与实际负载下的静态传递误差对比量级应一致第三层剥落尺寸趋于零时仿真结果应平滑过渡为健康工况的结果。这三层都过了程序才算真正可用。我刚写这套程序时也走了不少弯路最深刻的体会是剥落故障仿真成败的关键不在动力学求解器而在啮合刚度算得准不准。刚度曲线只要有5%的误差后续频谱边带特征就可能面目全非。建议任何刚入坑的朋友先从健康齿轮的刚度计算和验证做起把势能法公式和有限元结果对上了再往模型里加剥落参数这样定位问题会快很多。另外剥落的三个几何参数一定要定义成程序顶部的全局变量方便做参数敏感性扫描——当你需要画剥落长度从1 mm到5 mm变化时的响应瀑布图时就知道这个设计有多省事了。这组程序后续还可以扩展出齿根裂纹、齿面磨损等不同故障模型只需替换刚度计算模块中的故障几何描述动力学主循环完全不用动。