简介：这份MATLAB脚本面向齿轮动力学与机械振动方向的研究生、工程师及科研人员，用于求解直齿轮啮合过程中随时间动态变化的啮合刚度。核心算法基于势能法，通过建立轮齿变形能与外载荷之间的关系，结合齿面接触几何、材料弹性特性及啮合位置变化，逐点计算不同时刻的刚度数值，可生成单对齿啮入至啮出全过程的刚度曲线，输出时间序列对应的刚度数组，直接服务于齿轮系统动力学建模、振动响应分析与故障特征提取等后续研究。资源包共6个文件，约42KB，以m脚本为主体，辅以py脚本、txt说明、png结果图及工程配置文件，无需额外工具箱，在MATLAB基础平台即可运行。脚本输入参数涵盖模数、齿数、压力角、齿宽、弹性模量、泊松比及转速等基础齿轮参数，结构清晰，便于按需修改与二次开发。目前已有36人学习，适合希望快速获取时变啮合刚度计算工具、减少重复编程工作的读者参考使用。
1. 直齿轮啮合刚度为什么不能当常数：从一次振动超标排查说起
直齿轮啮合刚度随时间变化的MATLAB计算脚本，核心要解决的就是一件事：把齿轮副在啮合过程中刚度的周期性波动算出来。很多做齿轮箱振动分析的人一开始都会把啮合刚度当成一个定值代进动力学方程，结果仿真出来的频谱干干净净，实际测出来的振动信号却在啮合频率及其倍频附近出现一堆边带，怎么调阻尼都对不上。问题就出在刚度上——直齿轮重合度通常在1到2之间，单对齿啮合和双对齿啮合交替出现，参与承载的齿对数在变化，啮合刚度自然随时间周期性变化。这个时变特性正是齿轮系统参数振动的激励源，也是边频带的成因。
势能法实现的好处在于物理意义清晰：把轮齿当成变截面的悬臂梁，分别算弯曲势能、剪切势能、轴向压缩势能和接触势能，再对总势能求偏导得到刚度。相比有限元接触分析，势能法不需要建三维模型、不需要处理接触收敛，一个脚本几分钟就能跑完一个啮合周期，适合做参数扫描和动力学仿真的刚度输入。适合谁用？做齿轮动力学建模的研究生、做齿轮箱故障诊断的工程师、需要快速评估齿形参数对动态特性影响的设计人员。前提是你得接受它的假设：小变形、线弹性、齿根圆角简化处理，精度对常规模数齿轮够用，对修形齿轮或薄壁齿轮要谨慎。
2. 势能法算啮合刚度的力学模型与公式拆解
2.1 把轮齿拆成变截面悬臂梁：四种势能的物理含义
势能法的基本思路不复杂：轮齿受载后会发生弹性变形，变形过程中外力做功转化为弹性势能，刚度就是势能对变形的二阶导数。具体到直齿轮，一对啮合齿在接触点受力，力沿啮合线方向传递，轮齿可以等效为一根变截面悬臂梁，固定端在齿根圆处，自由端在齿顶。
在这个等效下，轮齿的弹性势能由四部分组成。弯曲势能来自弯矩作用，齿根附近截面小、弯矩大，贡献最大。剪切势能来自剪力作用，对短粗齿不可忽略。轴向压缩势能来自压力沿齿高方向的分量，通常占比最小但计算时不能直接扔掉。接触势能来自两齿面在接触点的局部变形，用赫兹接触理论处理，和齿面曲率半径、材料弹性模量、泊松比有关。
这四项势能分别对应四个刚度分量：弯曲刚度、剪切刚度、轴向压缩刚度、接触刚度。单对齿的综合刚度是这四个分量串联的结果，串联公式和电阻并联形式相反，是倒数相加再取倒数。这里有个容易翻车的地方：很多人直接把四个刚度相加，那就错了，串联是柔度相加。
注意：接触刚度在双齿啮合区要按两对齿的接触柔度并联处理，不是简单把单对齿刚度乘2。
2.2 单对齿刚度公式推导：从势能到刚度的关键一步
设啮合线方向的载荷为F，接触点沿啮合线方向的位移为δ，则刚度为k=F/δ。势能法不直接求位移，而是先写出势能U关于载荷F的表达式，然后利用卡氏定理δ=∂U/∂F，再代入刚度定义。
弯曲势能表达式为U_b=∫(M²/2EI)dx，其中M是弯矩，E是弹性模量，I是截面惯性矩。对直齿轮，截面是矩形，宽度为齿宽b，高度为2y（y是接触点到齿中心线的距离），所以I=b(2y)³/12。弯矩M=F·(h-x)，h是接触点高度，x是沿齿高方向的积分变量。积分从齿根圆角处积到接触点。
剪切势能U_s=∫(1.2F²/2GA)dx，1.2是矩形截面的剪切修正系数，G是剪切模量。轴向压缩势能U_a=∫(F²/2EA)dx。接触势能U_h=F²/(2k_h)，k_h是赫兹接触刚度，由接触半宽和材料参数决定。
把四项势能对F求偏导得到各分量柔度，再取倒数得到各分量刚度，最后串联。整个推导过程在脚本里体现为数值积分，因为截面高度y随x变化，解析积分对渐开线齿廓不现实。
2.3 重合度与啮合区划分：单齿区和双齿区怎么定
直齿轮重合度ε_α一般在1.2到1.8之间，意味着一个啮合周期里，有(2-ε_α)的比例是双齿啮合，(ε_α-1)的比例是单齿啮合。具体划分靠啮合线长度和基节的关系。
啮合线从主动轮齿根与从动轮齿顶接触点开始，到主动轮齿顶与从动轮齿根接触点结束。基节p_b=πm cosα，m是模数，α是压力角。啮合线长度L=√(r_a1²-r_b1²)+√(r_a2²-r_b2²)-a sinα，其中r_a是齿顶圆半径，r_b是基圆半径，a是中心距。
沿啮合线按基节划分，每过一个基节就有一对齿进入或退出啮合。脚本里用离散点遍历啮合线，每个点判断当前有几对齿在接触，分别算单对齿刚度再按并联处理。这里有个细节：双齿区两对齿的接触点位置不同，各自的弯曲力臂和接触曲率半径都不同，不能共用一套几何参数。
3. MATLAB脚本实现：从几何参数到刚度曲线的完整流程
3.1 齿轮基本参数定义与渐开线齿廓生成
脚本第一步是把齿轮参数写清楚，后面所有几何计算都依赖这些值。我一般把参数集中放在脚本开头，方便改模数、齿数、变位系数做扫描。
% 直齿轮基本参数 m = 3; % 模数 mm z1 = 20; % 主动轮齿数 z2 = 40; % 从动轮齿数 alpha = 20*pi/180; % 压力角 rad b = 20; % 齿宽 mm ha_star = 1; % 齿顶高系数 c_star = 0.25; % 顶隙系数 x1 = 0; % 主动轮变位系数 x2 = 0; % 从动轮变位系数 E = 206000; % 弹性模量 MPa nu = 0.3; % 泊松比 G = E/(2*(1+nu)); % 剪切模量 MPa % 分度圆、基圆、齿顶圆、齿根圆 d1 = m*z1; d2 = m*z2; db1 = d1*cos(alpha); db2 = d2*cos(alpha); da1 = d1 + 2*m*(ha_star + x1); da2 = d2 + 2*m*(ha_star + x2); df1 = d1 - 2*m*(ha_star + c_star - x1); df2 = d2 - 2*m*(ha_star + c_star - x2); a = (d1+d2)/2; % 标准中心距这段代码定义了齿轮的基本几何和材料参数。模数和齿数决定尺寸，压力角影响基圆，变位系数改变齿顶齿根圆。弹性模量和泊松比用于算接触刚度，剪切模量用于剪切势能。注意单位统一用毫米和兆帕，后面积分时长度单位是毫米，力单位是牛，刚度单位是牛每微米，这是齿轮动力学里常用的量纲。
渐开线齿廓生成用参数方程，从基圆到齿顶圆取一系列展开角。脚本里不需要把齿廓画出来，但需要计算任意接触点对应的齿厚和截面高度，所以齿廓的几何关系要写对。
% 渐开线齿廓参数：从基圆展开到齿顶圆 r_k = linspace(db1/2, da1/2, 200); % 主动轮齿廓半径序列 alpha_k = acos(db1./(2*r_k)); % 各点压力角 inv_alpha = tan(alpha_k) - alpha_k; % 渐开线函数 % 齿厚半角 s_half = (pi/(2*z1) + 2*x1*tan(alpha)/z1) + (tan(alpha) - alpha) - inv_alpha; y_k = r_k .* sin(s_half); % 接触点距齿中心线的距离 h_k = r_k .* cos(s_half) - db1/2; % 接触点距基圆的径向高度这里y_k是接触点到齿中心线的垂直距离，对应悬臂梁截面半高。h_k是接触点相对基圆的径向高度，用于算弯矩力臂。齿厚半角公式里包含了变位项，标准齿轮x1=0时退化为常规公式。200个离散点够用，再多对精度提升有限但计算量线性增长。
3.2 四种势能的数值积分实现
有了齿廓几何，接下来对每个接触点位置算四项势能。积分变量是沿齿高方向的坐标，从齿根过渡曲线终点积到接触点。齿根过渡曲线用圆弧近似，起点取齿根圆与渐开线交点。
function [Ub, Us, Ua] = tooth_energy(F, r_k, y_k, h_k, b, E, G) % 输入：F啮合力，r_k齿廓半径序列，y_k截面半高，h_k径向高度 % 输出：弯曲、剪切、轴向压缩势能 n = length(r_k); Ub = 0; Us = 0; Ua = 0; for i = 1:n-1 dx = r_k(i+1) - r_k(i); % 积分步长 y_mid = (y_k(i)+y_k(i+1))/2; % 中点截面半高 h_mid = (h_k(i)+h_k(i+1))/2; % 中点力臂 I = b*(2*y_mid)^3/12; % 截面惯性矩 A = b*2*y_mid; % 截面积 M = F*h_mid; % 弯矩 Ub = Ub + M^2/(2*E*I)*dx; Us = Us + 1.2*F^2/(2*G*A)*dx; Ua = Ua + F^2/(2*E*A)*dx; end end这个函数对每个微元段用中点矩形法积分。弯曲势能里弯矩M=F*h_mid，力臂是接触点相对当前截面的径向高度差，这里简化用h_mid近似，严格应该用h_k(i)-h_k(j)的差值，但齿根到接触点的径向距离远大于截面高度变化，简化误差在5%以内。剪切修正系数1.2是矩形截面标准值，圆形截面是10/9，齿轮齿截面接近矩形所以用1.2。
接触势能单独算，用赫兹线接触公式。接触刚度k_h与接触点曲率半径和材料参数有关。
function kh = contact_stiffness(F, rho1, rho2, b, E, nu) % 赫兹接触刚度 E_star = E/(2*(1-nu^2)); % 等效弹性模量 R_eq = 1/(1/rho1 + 1/rho2); % 等效曲率半径 % 接触半宽 a_h = sqrt(4*F*R_eq/(pi*b*E_star)); % 接触刚度（线接触） kh = pi*b*E_star/(4*(log(4*R_eq/a_h)-0.5)); endrho1和rho2分别是两齿在接触点的曲率半径，由渐开线性质算得。接触刚度随载荷变化，因为接触半宽a_h与F的平方根成正比，所以整个啮合刚度曲线其实是非线性的。做线性动力学分析时通常在额定载荷附近取切线刚度，做非线性分析才保留这个关系。
3.3 单双齿交替区的刚度合成与曲线绘制
单对齿综合刚度是四个分量串联，双齿区是两对齿的柔度并联。脚本主循环遍历啮合线，每个位置判断接触齿对数。
% 啮合线离散 N = 500; s = linspace(0, L_mesh, N); % L_mesh为啮合线长度 k_total = zeros(1, N); for i = 1:N % 判断当前接触齿对数 n_pairs = count_pairs(s(i), p_b, L_mesh); k_pair = zeros(1, n_pairs); for j = 1:n_pairs % 算第j对齿的接触点几何 [rho1, rho2, h_contact] = contact_geometry(s(i), j, ...); % 算四项势能对应的柔度 [Ub, Us, Ua] = tooth_energy(1, ...); % 单位力 delta_b = 2*Ub; delta_s = 2*Us; delta_a = 2*Ua; kh = contact_stiffness(1, rho1, rho2, b, E, nu); delta_h = 1/kh; delta_total = delta_b + delta_s + delta_a + delta_h; k_pair(j) = 1/delta_total; end % 并联：柔度相加再取倒数 k_total(i) = 1/sum(1./k_pair); end plot(s, k_total, 'LineWidth', 1.5); xlabel('啮合线位置 (mm)'); ylabel('啮合刚度 (N/\mum)');count_pairs函数根据当前啮合线位置和基节判断有几对齿在接触。接触点几何函数算该对齿在当前啮合位置的曲率半径和力臂。注意势能函数传单位力1，因为柔度就是单位力下的位移，这样算完直接得柔度，取倒数就是刚度。
跑出来的曲线典型特征是：双齿区刚度高且平缓，单齿区刚度低且呈凹形，过渡点有阶跃。阶跃处是齿对进入或退出啮合的瞬间，实际由于齿顶修缘和制造误差会平滑一些，但理想模型就是阶跃。这个阶跃正是动力学仿真里参数振动的激励。
4. 避坑与排查：势能法算啮合刚度最容易翻车的五个地方
4.1 刚度曲线出现负值或突变尖峰
现象：跑出来的k_total在某些点变成负数，或者出现不合理的尖峰。原因通常是积分区间搞反了，齿根到接触点的积分方向如果从接触点积到齿根，dx为负，势能累加就出错。另一个可能是接触点位置超出了齿廓有效范围，比如算到了基圆以下，那里没有渐开线齿廓，y_k算出虚数或零，截面惯性矩为零导致柔度无穷大。
解决：在积分前加判断，接触点半径必须大于基圆半径加一个安全余量。积分方向统一从齿根到接触点，dx取绝对值。如果接触点接近齿根，弯曲力臂趋近零，弯曲柔度也趋近零，这是合理的，但数值上要防止除零。
4.2 双齿区刚度反而比单齿区低
现象：双齿啮合区算出来的综合刚度比单齿区还小，和物理直觉相反。原因多半是并联公式写错了，把柔度并联写成了刚度并联。双齿区两对齿同时承载，总柔度是两对齿柔度并联，总柔度变小，总刚度变大。如果代码里写成k_total=k1+k2再取倒数，那就反了。
解决：记住口诀“刚度串联取倒数和，柔度并联直接加”。单对齿内部四个分量是串联，用柔度相加；双齿区两对齿是并联，也是柔度相加。整个链条里只有最后一步取倒数得到刚度。
4.3 重合度算错导致单双齿区比例不对
现象：刚度曲线里单齿区和双齿区的宽度比例明显不对，比如单齿区占了80%。原因通常是重合度公式用错，或者啮合线长度算错。重合度ε_α=L_mesh/p_b，如果L_mesh少算了一个齿顶高，ε_α偏小，单齿区就偏大。
解决：用标准公式验算重合度，ε_α=[z1(tanα_a1-tanα)+z2(tanα_a2-tanα)]/(2π)，其中α_a是齿顶压力角。两种方法算出来应该一致。啮合线起点和终点的判断要仔细，起点是主动轮齿根和从动轮齿顶接触，终点反过来。
4.4 接触刚度数量级异常
现象：接触刚度算出来比弯曲刚度小两个数量级，导致总刚度几乎全由接触刚度决定。原因通常是接触半宽公式里单位混了，比如力用牛但长度用米，或者弹性模量用兆帕但长度用毫米，量纲不一致。
解决：全程用毫米-牛-兆帕体系。接触半宽a_h的单位是毫米，等效曲率半径R_eq单位毫米，线载荷F/b单位牛每毫米。算完检查一下a_h是否在0.01到0.1毫米量级，接触刚度是否在10^5到10^6牛每毫米量级，偏离太多就是单位问题。
4.5 数值积分步长导致曲线毛刺
现象：刚度曲线不平滑，有高频毛刺，尤其单齿区凹面不光滑。原因是齿廓离散点太少，或者积分步长太大，中点矩形法精度不够。200个点对模数3的齿轮够用，模数小或者齿数多的时候齿廓曲率变化快，需要加密。
解决：把齿廓离散点增加到500以上，或者改用梯形积分。更稳妥的做法是自适应步长，在齿根附近加密，齿顶附近稀疏。跑完看曲线一阶导是否连续，毛刺严重就加密重跑。
5. 从刚度曲线到动力学仿真：时变刚度的调用与验证技巧
拿到啮合刚度曲线只是第一步，真正要用起来得把它喂给动力学方程。直齿轮扭转振动模型里，啮合刚度是时变系数，方程是Mθ''+Cθ'+K(t)θ=F。K(t)就是刚算出来的k_total，按啮合位置映射到时间域。映射关系是啮合线位置s=ω·r_b·t，ω是角速度，r_b是基圆半径。这里有个细节：主动轮和从动轮的基圆半径不同，但啮合线速度相同，用哪个算都一样。
我一般把刚度曲线先做傅里叶级数展开，取前5到8阶谐波，因为动力学仿真里高阶谐波对响应影响小，但保留太多阶会让步长被迫取很小。展开后K(t)=k_mean+Σk_i·cos(i·ω_m·t+φ_i)，ω_m是啮合频率。均值k_mean就是平均啮合刚度，各阶谐波幅值决定边频带强度。
验证方法有两个。一是能量守恒检查：在一个啮合周期内，刚度曲线积分得到的平均刚度应该和用ISO标准公式算的啮合刚度接近，偏差超过15%说明模型有问题。二是和有限元结果对比：如果手头有接触分析结果，取相同载荷和边界条件，刚度曲线形状应该一致，峰值位置对应单双齿交替点。
一个实用技巧：做参数扫描时，把齿廓离散点数和啮合线离散点数做成可配置参数，先粗跑看趋势，找到敏感区间再细跑。我习惯先用100个齿廓点加200个啮合线点跑一遍，看刚度均值和谐波幅值，如果和预期差得远，先查几何参数而不是急着加密网格。另一个习惯是每次改完参数把重合度打印出来，ε_α低于1.1或高于1.9都要警惕，前者接近单齿啮合极限，后者接近三齿啮合，势能法的两齿假设就不成立了。
最后说个血泪教训：势能法算的刚度是静刚度，动力学仿真里如果转速很高，还得考虑动态刚度的频变特性，但那是另一个话题了。先把静刚度算准，曲线形状对了，边频带位置就不会错。希望帮到你。
本文还有配套的精品资源，点击获取
