企业官网建设流程全解析

做齿轮动力学或者故障诊断的朋友，大概率都绕不开时变啮合刚度（TVMS）这个参数。很多论文里都会画一条波浪形的刚度曲线，说这是齿轮副的“内部激励源”。但真到自己动手算的时候，尤其是想把齿根裂纹这种局部故障也塞进模型里，就会发现好多细节没写清楚。我之前就是照着马辉、罗阳等文献的思路，用MATLAB写了一个考虑齿根裂纹的直齿轮时变啮合刚度计算程序，中间踩了不少坑，也反复对比过有限元结果。这篇博文就把整个计算思路、代码架构、裂纹建模逻辑，以及那些文献里不太会写明的实操细节，系统整理出来。

这套方法的核心是势能法，思路是把轮齿当成变截面悬臂梁，把啮合力引起的能量分成弯曲、剪切、轴向压缩、赫兹接触和基体弹性变形几个部分，再通过能量守恒反推刚度。齿根裂纹的影响，最终落到截面积和截面惯性矩的折减上。整个过程不涉及任何商业软件，用MATLAB矩阵化编程完全可以跑，非常适合做参数敏感性分析，比如裂纹深度从10%到80%逐级变化时，刚度曲线到底怎么变、振动响应会有什么趋势，这类工作用有限元挨个建模会非常痛苦，而势能法编译一次就能批量出结果。

1. 为什么选势能法，以及时变啮合刚度的本质

1.1 时变啮合刚度到底在算什么

先明确一件事：齿轮啮合过程中，接触点是沿啮合线移动的，轮齿的受力位置和力臂时刻在变。再加上重合度通常不是整数，所以啮合过程会出现单齿啮合区和双齿啮合区交替的规律。双齿区里两对齿同时分担载荷，总变形量反而比单齿区小，刚度也就更高。这个随啮合时间或转角变化的综合刚度，就是时变啮合刚度。

从工程意义上看，这个刚度波动是齿轮系统最核心的动态激励源之一，它直接决定了齿轮振动的幅值和频率成分。齿根出现裂纹后，局部柔性增大，刚度曲线会在裂纹齿参与啮合的位置出现一个明显的凹陷，这个凹陷正是故障诊断里“边频调制”的机理来源。所以算准TVMS，不管是做动力学响应预报、振动信号仿真，还是做损伤识别特征提取，都是第一步。

1.2 势能法相比有限元法的核心优势

最初考虑过用有限元软件扫参数，后来放弃了。原因很简单：裂纹参数一改，就得重新建模、重新剖分网格，尤其是裂纹尖端的网格要加密，几千个样本跑下来，光网格处理的时间就让人崩溃。而势能法把轮齿抽象成悬臂梁，裂纹的影响只是改变了积分里的截面参数，算一次基准模型只需要毫秒级，批量扫描深度、角度、齿数等参数非常合适。

当然，势能法也有局限。它本质上是二维平面的简化模型，对齿面接触的局部弹性问题处理得比较粗糙，应力集中系数也需要经验修正。所以我的态度是：做了几百组参数对比和趋势分析，然后挑几个典型工况用二维有限元做交叉验证，两边印证着来，既保留效率，又不至于偏离真实物理太远。

1.3 从马辉、罗阳等文献中提炼的建模逻辑

马辉、罗阳等文献在含裂纹齿轮刚度计算上的处理方式，基本可以概括成三步：

1. 把单个轮齿等效成固定于齿根圆的变截面悬臂梁，齿廓渐开线部分离散成若干截面；
2. 对每一个微小截面计算面积、惯性矩等几何参数；
3. 将齿根裂纹视为截面有效面积的削减，把含裂纹的几何参数重新代入弯曲、剪切刚度的积分表达式。

这个过程简洁但严谨。文献里对裂纹简化的前提通常是贯穿齿宽的直线裂纹，也就是把三维裂纹问题先压成二维来处理。这篇博文的代码也是这样假设的，因为工程上高周疲劳的齿根裂纹早期往往呈线状扩展，先贯穿齿宽再向内扩展的简化是合理的。

2. 齿根裂纹如何进入刚度计算模型

2.1 裂纹几何参数怎么定义

裂纹模型需要参数化，我用的两个核心参数是裂纹深度和裂纹角度。

深度我用相对值表示，比如q表示裂纹尖端沿垂直于齿体表面方向扩展的深度占全齿高的比例。角度v表示裂纹扩展方向与齿体中心线的夹角，通常情况下裂纹从齿根圆角应力最大点出发，沿与齿面法线呈一定角度的方向扩展。文献里为了计算方便，普遍假设裂纹是一条直线，这样每个截面上被削掉的部分就是线性变化的。

实际使用中，早期故障推荐从q=5%~10%起步，因为在这个范围内刚度下降很小，这正好对应“早期裂纹难以诊断”这一工程现象。深度超过30%之后，刚度变化就会变得非常明显，振动特征也开始突出。

2.2 截面惯性矩的折减原理

这是整个裂纹建模的核心。

轮齿的渐开线齿廓在不同高度处的齿厚不同，悬臂梁模型里每个截面x位置都有自己的厚度S_x。一旦出现裂纹，裂纹尖端以下的材料仍然起支撑作用，但裂纹尖端以上部分会形成“断开”区域，原本有效的抗弯截面被削弱。计算时，我直接改写每个x截面上的有效厚度S'_x，如果裂纹把该截面分成了两段，则只保留仍与齿体主体相连的那一段。

更准确的做法是：先由裂纹直线方程求出它与齿体轮廓的交点，判断该截面上裂纹覆盖的宽度范围，再用数值方法求出剩余截面面积和惯性矩。注意，这里不能用简单的“厚度乘以一个系数”代替，因为惯性矩和厚度的三次方成正比，同样的厚度削减比例，惯性矩损失更大，对刚度的影响也显著得多。

2.3 不同刚度分量对裂纹的敏感程度

我实际算下来，不同能量分量对裂纹的敏感度差异很大，这直接影响故障特征的解释。

• 弯曲刚度：最敏感。因为弯曲变形正比于力臂的贡献和截面惯性矩的倒数，裂纹导致惯性矩下降，弯曲柔度显著增加。
• 剪切刚度：比较敏感。剪切变形取决于截面积，裂纹使有效截面积下降，刚度也会降低，但影响幅度通常比弯曲小。
• 轴向压缩刚度：基本可以忽略。齿面法向力分解出的轴向分量较小，且轴向压缩刚度由整个截面面积决定，裂纹对面积的影响有限。
• 赫兹接触刚度：不随裂纹变化。它只由接触点的曲率半径和材料弹性模量决定，不涉及裂纹几何。
• 基体弹性刚度：几乎不变。这个分量反映轮体基体对齿根支承的柔度，裂纹在齿体局部，对基体整体的弹性影响很小。

这个结论很重要，后续做故障诊断特征提取时，主要盯弯曲刚度的缺口来判断裂纹程度，而不是笼统看总刚度。

3. MATLAB实现详解：从参数定义到刚度曲线输出

3.1 计算流程总览

整个程序的执行顺序推荐这样安排：

1. 输入齿轮基本参数（模数、齿数、压力角、齿宽、材料参数）；
2. 计算齿廓几何，确定基圆、分度圆、齿顶圆、齿根圆以及渐开线离散点；
3. 计算重合度，确定一个啮合周期内单双齿啮合区间的边界；
4. 对啮合线上每一个离散位置，计算该啮合点相对齿根的位置和力臂；
5. 分别计算无裂纹和含裂纹两种情况下的五种刚度分量；
6. 按单双齿状态组合成总时变啮合刚度；
7. 绘图输出，并对关键工况做验证。

建议把五种刚度的计算封装成独立子函数，便于替换和调试。另外程序全程保持单位一致，几何量用mm，力用N，弹性模量用MPa（即N/mm²），算出来的刚度单位就是N/mm。

3.2 齿轮几何参数计算

先给定一个标准算例的参数，方便后续对照：

% 基本参数 m = 2; % 模数 mm z1 = 20; % 小齿轮齿数 z2 = 30; % 大齿轮齿数 alpha = 20*pi/180; % 压力角 rad B = 20; % 齿宽 mm E = 2.06e5; % 弹性模量 MPa (N/mm^2) nu = 0.3; % 泊松比 r1_p = m*z1/2; % 分度圆半径 r2_p = m*z2/2; r_b1 = r1_p*cos(alpha); % 基圆半径 r_b2 = r2_p*cos(alpha); r1_a = r1_p + m; % 齿顶圆半径 r2_a = r2_p + m; r_f1 = r1_p - 1.25*m; % 齿根圆半径 r_f2 = r2_p - 1.25*m;

这里有个新手容易犯的错误：在计算变截面悬臂梁积分下限时，有些程序直接从基圆算起，这是不对的。有效的啮合起始点应该是配对齿轮齿顶圆与该齿轮渐开线相交的点，而不是基圆本身。如果直接从基圆作为积分起点，算出来的刚度会偏小，导致曲线整体偏低。正确做法是先计算理论啮合线长度和啮合起始点半径。

重合度用下式计算：

g_a = sqrt(r1_a^2 - r_b1^2) + sqrt(r2_a^2 - r_b2^2) - (r1_p + r2_p)*sin(alpha); epsilon = g_a / (pi*m*cos(alpha));

这个重合度决定了单双齿分界。比如epsilon=1.6，就说明啮合周期内约60%是双齿区，40%是单齿区，具体边界点要换算成主动轮的转角。

3.3 五种刚度分量的计算子程序

这里把每个子程序的关键逻辑说一下，代码格式可以直接照着用。

赫兹接触刚度是常数，与位置无关。平面应变状态下公式为：

function kh = hertz_stiffness(E, B, nu) kh = pi*E*B / (4*(1 - nu^2)); end

这个公式里的系数取决于接触模型，对于两个齿廓的线接触，按无限长圆柱体赫兹接触导出，系数通常是pi/4。注意平面应力状态下系数会不一样，齿轮齿宽足够大时按平面应变处理更符合实际。

弯曲刚度是整个程序里最核心的部分。它的柔度表达式是：

1/k_b = ∫ 0^d [F_b*(d-x) - F_ah(x)]^2 / (EI(x)) dx

其中d是啮合点沿齿高方向到齿根的距离，h(x)是截面到力作用线的偏移量，I(x)是截面惯性矩。含裂纹时，I(x)要换成I'(x)。

function kb = bending_stiffness(params, x_mesh, I_x, F_b, F_a) d = params.d_mesh; h_x = params.h_mesh; integrand = (F_b*(d - x_mesh) - F_a*h_x).^2 ./ (params.E .* I_x); kb = 1 / trapz(x_mesh, integrand); end

积分用trapz做数值积分就可以了，不需要上符号积分。关键是x_mesh的离散密度，我建议一个啮合周期至少取300个位置点，每个齿轮廓在积分方向的离散点也要在100以上，太疏了曲线会有锯齿。

剪切刚度的表达式为：

1/k_s = ∫ 0^d 1.2F_a^2 / (GA(x)) dx

这里系数1.2是矩形截面的剪切修正系数。G = E / (2*(1+nu))。同样，含裂纹时A(x)要替换成有效截面积。

轴向压缩刚度公式：

1/k_a = ∫ 0^d F_a^2 / (E*A(x)) dx

实际计算时这个分量占比很小，但仍建议保留，因为严谨性需要。F_a是法向力分解出的轴向分量。

基体弹性刚度我采用Sainsot等文献给出的拟合公式，处理轮体基体柔度：

1/k_f = cos^2(alpha) / (EB) * [ L(u_f/S_f)^2 + M*(u_f/S_f) + P*(1 + Q*tan^2(alpha)) ]

这里L、M、P、Q是拟合系数，u_f和S_f由齿根圆半径、齿根圆弧和齿厚计算得到。这个公式看着繁琐，但好处是不用建轮体模型就能考虑基体柔度，效率很高。编写时为这几个几何量单独写一个计算函数会清爽很多。

3.4 时变啮合刚度主循环与单双齿组合

将齿轮副一个啮合周期内各啮合位置的刚度拼装起来，是整个程序的主线。这里有个细节：啮合力作用线和齿廓渐开线的角度在不同位置略有不同，但很多程序为了简化，直接按标准压力角处理。对于精度要求不高的情况可以接受，如果要更严谨，就应该在循环内重新计算瞬时啮合角。

单齿区和双齿区的组合方式：齿对1和齿对2是并联关系，总弹性变形等于各对齿弹性变形之和。所以双齿区总柔度等于两对齿柔度相加，再取倒数得到总刚度：

k_total = 1 / (1/k_pair1 + 1/k_pair2)

而单齿区直接取那一对齿的总刚度即可。

主循环的大致逻辑：

% 归一化转角：0到1为一个啮合周期 N = 501; % 啮合位置采样数 theta_mesh = linspace(0, 2*pi/ z1, N); k_total = zeros(1, N); for i = 1:N % 根据转角确定当前啮合位置，计算齿对1的啮合点半径 r_mesh1 = calc_mesh_radius(theta_mesh(i)); % 判断是否处于双齿区 if in_double_contact(theta_mesh(i), epsilon) % 齿对1 + 齿对2 k_pair1 = calc_pair_stiffness(r_mesh1, crack_params, ...); k_pair2 = calc_pair_stiffness(r_mesh2, crack_params, ...); k_total(i) = 1 / (1/k_pair1 + 1/k_pair2); else k_pair1 = calc_pair_stiffness(r_mesh1, crack_params, ...); k_total(i) = k_pair1; end end

绘图部分我习惯把无裂纹和有裂纹的曲线画在同一张图里对比，这样故障的影响范围一目了然：

figure('Color','w'); plot(theta_mesh*180/pi, k_health, 'k-', 'LineWidth', 1.5); hold on; plot(theta_mesh*180/pi, k_crack, 'r--', 'LineWidth', 1.5); xlabel('主动轮转角 (deg)'); ylabel('时变啮合刚度 (N/mm)'); legend('无裂纹', '含裂纹'); grid on;

从我实际跑出来的结果看，裂纹深度30%时，刚度曲线在裂纹齿对应区间会出现约8%~15%的局部下降，下降幅度和齿数、重合度有关。如果裂纹放在从动轮上，凹坑的位置会偏移，这个相位信息在故障定位里非常有用。

4. 常见问题排查、验证方法与实操避坑

4.1 刚度曲线形状异常先查几何边界

如果出来的曲线单双齿过渡点位置不对，或者双齿区的刚度比单齿区还低，几乎可以断定是啮合区间划分或者啮合起始点出问题了。排查顺序：

1. 先核对重合度计算结果。手算一遍，确认程序里sqrt(r_a^2 - r_b^2)这部分取的是同一齿轮的数值，别把主动轮和从动轮的半径混用。
2. 再核对啮合起始点半径。正确表达式是由配对齿轮齿顶圆确定的，不是基圆。
3. 最后检查单双齿区间的边界点换算。要把啮合线长度比例换算成主动轮的转角范围，这个换算用到基圆半径，别用分度圆半径。

我第一版程序就是栽在第三个问题上，出来的双齿区宽度明显不对，反复排查才发现是转角换算写错了。

4.2 结果可信度怎么验证

没有裂纹的模型最好验证。找几篇经典文献里的标准齿轮算例，把参数输进去，对比平均啮合刚度和刚度波动幅度。一般来说，解析法和势能法结果差异在5%以内是正常的。如果差异大，优先检查基体柔度项是否没加，以及齿根支撑位置是否取错了。

含裂纹的结果验证稍微麻烦些。可以拿有限元二维模型做几个典型深度点对标，比如10%、30%、50%三个深度，对比刚度下降百分比。我的经验是趋势一致就算合格，绝对值的误差控制在10%以内已经很理想了。当裂纹深度接近穿透齿根时，势能法计算会出现数值不稳定，因为悬臂梁假设在极端损伤下已经失真，这个区间就不要强行用它了。

4.3 几个容易忽略但影响很大的细节

第一个细节是齿根过渡圆角。很多初始版本程序直接把齿根支撑点放在齿根圆上，完全不考虑过渡圆角，这样算出来的基体刚度会偏大。处理办法是把支撑点沿过渡圆弧稍微向内移动一个距离，或者用等效法把过渡圆角的影响折进齿根厚度里。

第二个细节是矩阵化加速。单个子程序用for循环慢慢算也能跑，但参数扫描的时候要跑几百次，速度差别就出来了。建议对齿廓离散、积分计算做向量化处理，尽量用trapz替代复杂的循环累加。

第三个细节是裂纹表达式的连续性。裂纹深度从0逐渐增加时，刚度曲线应当平滑过渡。如果发现曲线明显跳变，八成是裂纹与齿廓相交的判断条件写得不连续，检查逻辑分支有没有覆盖所有几何可能性。

第四个细节是绘图线型。我习惯无裂纹用实线看全局，裂纹用虚线叠加，对比时能明显看到凹陷位置和深度。线宽建议设置1.5以上，否则导出到论文里会显得很淡。

4.4 这个模型的扩展场景

算出来的时变啮合刚度可以直接接到MATLAB的动力学求解流程里，用ode45求解齿轮副的扭转振动模型，就能得到考虑故障的动态响应。再进一步，把刚度曲线的下降幅度作为故障特征，可以做不同裂纹深度的模式识别。

另一个方向是把刚度结果用于裂纹扩展寿命预测。虽然势能法不能直接算出裂纹尖端的应力强度因子，但可以从刚度退化反推载荷幅度，结合Paris公式估算剩余寿命，这在状态检修里很有工程价值。

如果后面要往斜齿轮上扩展，思路也不变，只是把接触线从“一条直线”变成“斜线”，需要把齿宽方向分成多个薄片，每个薄片按直齿轮处理，再叠加求总刚度。代码架构上，只需要把主循环改成双重循环，外层遍历齿宽切片，内层走原来的单齿计算逻辑。

写在最后的实操心得

这个程序前前后后我改过三版，最大的体会是：势能法本身公式并不复杂，难点全在齿轮几何的边界条件上。尤其含裂纹的时候，每个截面的有效面积和惯性矩都要仔细判断，稍不注意就会出现不连续点。代码里每个几何量我都建议把公式来源写到注释里，不然隔几周回来看就容易懵。

另外，算出来的刚度曲线不要只看总刚度，最好把弯曲、剪切、基体这几项分别画出来观察。裂纹对弯曲项的削弱最明显，如果总曲线变化不大，先看弯曲项有没有真的下降，这能帮你快速定位是模型问题还是裂纹参数设置问题。

你要是也在做齿轮故障诊断或者动力学仿真，建议把这个程序当作一个基础工具，先把无裂纹模型校准准确，再逐步加入裂纹参数。希望这些经验能帮你少走一些弯路。