企业官网建设流程全解析

搞齿轮故障诊断这些年，经常被同行问到的一句话是：剥落故障在振动信号上到底长什么样？说实话，光靠现场测振谱很难讲清楚，因为传感器采到的响应是齿轮副时变啮合刚度、齿侧间隙、误差激励、负载波动等多因素耦合后的综合结果。我一直建议做故障诊断的工程师把“含剥落故障直齿轮啮合刚度”的计算先写成程序，再把刚度结果喂进齿轮非线性动力学方程里做数值仿真，这样才能把故障特征从一堆耦合信号里单独拧出来看清楚。本文就从啮合刚度的机理推导开始，讲到非线性动力学模型的搭建和MATLAB程序实现，最后附上我实际调试中踩过的坑，希望能给正在做齿轮动力学仿真的朋友一些参考。

这篇文章适合三类人：一是做齿轮箱故障诊断、想在仿真里复现剥落特征的研究生和工程师；二是做传动系统设计、需要评估损伤对动态性能影响的机械工程师；三是刚接触齿轮非线性动力学、想把论文里的方程变成能跑出结果的程序的初学者。我会尽量把公式背后的物理含义讲透，代码结构也会拆开说明。

1. 剥落故障对啮合刚度的作用机理

1.1 剥落不是“表面刮花”，它直接改写齿轮副的承载路径

齿轮剥落（spalling）属于接触疲劳失效，本质是齿面次表层在循环接触应力作用下萌生裂纹，裂纹扩展到一定深度后引起材料成片脱落，形成不规则凹坑。很多人对剥落的理解停留在“齿面多了个坑”，但从动力学仿真的角度看，这个坑的存在意味着两个重要变化：第一，齿轮副在啮合到剥落区域时，原本由整个齿宽承担的接触载荷被迫转移，实际有效接触面积减小；第二，轮齿悬臂梁模型的截面几何发生变化，弯曲刚度、剪切刚度随之改变。这两点最终都归结到一个核心物理量——啮合刚度的时变曲线出现了局部的凹坑和突变。

理解这一点是后续一切工作的前提。齿轮啮合刚度并不是恒定值，它随啮合位置周期变化，这本身就是齿轮系统最主要的内部激励源。当剥落故障存在时，相当于在原本周期变化的刚度曲线上叠加了一个“局部缺陷脉冲”，这个脉冲会激励起齿轮系统的共振响应，在频谱上表现为啮合频率及其谐波附近的边带。我见过不少初学者直接拿恒定刚度建模，仿真出来的频谱干干净净，根本复现不出边带特征，就是因为没有把剥落引起的刚度变化放进去。

1.2 恒定刚度假设在故障仿真里的局限

经典齿轮设计中常用平均啮合刚度或ISO标准推荐的刚度公式来算承载能力，这在强度校核阶段够用，但做故障诊断仿真就远远不够了。齿轮副内部激励的主要来源就是时变啮合刚度（TVMS, Time-Varying Mesh Stiffness），恒定刚度等于把最大的激励源直接抹掉了。剥落故障诊断的核心逻辑是“故障改变了刚度→刚度改变激励→激励改变响应”，如果第一步就用了恒定刚度，后面全是空中楼阁。

1.3 程序化解决的核心思路

把剥落故障纳入动力学仿真的完整链路是：建立剥落的几何参数模型 → 基于势能法计算含故障的时变啮合刚度 → 把刚度序列作为时变系数代入非线性运动微分方程 → 数值积分得到振动响应 → 对响应做频谱、相图、庞加莱截面分析。这套链路每一步都有数学依据，每一步都可以用程序固化下来。我下面按这个顺序逐层展开。

2. 剥落故障的数学化描述与几何建模

2.1 矩形剥落模型的三要素

工程中为了解析计算方便，通常把不规则剥落坑简化为矩形凹槽。这个简化在故障诊断学术界是通用做法，McGraw-Hill的教材和大量论文里都用这种近似。矩形剥落需要定义三个参数：剥落轴向长度 (L)、剥落深度 (h_s)、剥落沿齿廓方向的位置（通常用距齿顶的距离或对应滚动角描述）。

我以自己常用的参数为例：模数 (m=3)，小齿轮齿数 (z_1=20)，大齿轮齿数 (z_2=30)，齿宽 (b=20) mm，压力角 (\alpha=20^\circ)，标准安装距。剥落参数取轴向长度 (L=3) mm、深度 (h_s=0.5) mm，位置在小齿轮单双齿啮合交替区域附近。剥落深度不能超过齿根危险截面的合理范围，否则刚度会算成负值，这在第7章会专门讲。

2.2 啮合过程的接触点移动轨迹

要计算时变啮合刚度，先得弄清楚一个啮合周期内接触点的位置怎么变化。直齿轮副啮合时，接触线沿齿廓方向从齿根向齿顶（或反过来）扫过，同时啮合的齿对数在1和2之间切换。啮合刚度的周期就是端面基节 (p_{bt} = \pi m \cos\alpha) 对应的小齿轮转角 (\varphi_z = 2\pi / z_1)。

程序实现时，我把小齿轮转过一个齿距角均分成 (N=200) 个离散位置，每个位置分别计算当前参与啮合的齿对数、每个齿上接触点的坐标，再根据接触点是否落入剥落区间来决定是否在刚度公式中计入剥落影响几何参数。这里有个关键细节：剥落位置若发生在单齿啮合区，刚度下降会非常剧烈，因为那时只有一对齿承担全部载荷；若在双齿啮合区，载荷由两对齿分担，刚度下降相对温和。这个位置敏感性本身就是故障诊断里区分剥落严重程度的重要依据。

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])（(i=1,...,N)，(N=400)），通过线性插值得到任意时刻的刚度值。误差激励 (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，我一般设RelTol=1e-6, AbsTol=1e-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) 处出现明显的边带簇，(k=1,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_g=0)）的条件下，刚度曲线和固有频率应与解析解或有限元结果误差在5%以内；第二层，把仿真位移响应的平均值与实际负载下的静态传递误差对比，量级应一致；第三层，剥落尺寸趋于零时，仿真结果应平滑过渡为健康工况的结果。这三层都过了，程序才算真正可用。

我刚写这套程序时也走了不少弯路，最深刻的体会是：剥落故障仿真成败的关键不在动力学求解器，而在啮合刚度算得准不准。刚度曲线只要有5%的误差，后续频谱边带特征就可能面目全非。建议任何刚入坑的朋友，先从健康齿轮的刚度计算和验证做起，把势能法公式和有限元结果对上了，再往模型里加剥落参数，这样定位问题会快很多。另外，剥落的三个几何参数一定要定义成程序顶部的全局变量，方便做参数敏感性扫描——当你需要画剥落长度从1 mm到5 mm变化时的响应瀑布图时，就知道这个设计有多省事了。这组程序后续还可以扩展出齿根裂纹、齿面磨损等不同故障模型，只需替换刚度计算模块中的故障几何描述，动力学主循环完全不用动。