做齿轮动力学仿真这几年我最大的感受是很多问题一开始并不是动力学问题而是刚度没算对。前阵子帮一家减速机厂家排查异响拆开一看齿面接触斑点明显偏向齿顶测振数据里啮合频率及其边频带着实扎眼。要解释清楚这个现象绕不开的就是时变啮合刚度。它随啮合位置周期性变化是齿轮系统最核心的激励源。而工程上算这个量我优先用的还是石川算法——模型简洁、物理意义清楚、MATLAB实现起来不费劲几分钟就能跑完一个参数工况比动辄几万节点的有限元模型实用得多。这篇文章就把整套实战过程写出来从石川算法的力学原理到MATLAB里的几何预处理、刚度积分、齿面接触变形和载荷分担再到把刚度曲线喂进单自由度动力学模型看振动响应最后聊聊故障齿轮的刚度特征和频谱识别。适合正在做齿轮动力学课题的学生、减速机设计工程师以及想用MATLAB快速评估齿轮修形或故障影响的同行参考。1. 为什么齿轮动力学分析绕不开时变啮合刚度1.1 刚度激励才是齿轮振动的源头齿轮传动中轮齿承受法向载荷后必然产生弹性变形。这个变形量由轮齿自身的几何形状和材料属性决定换算成力与变形的关系就是啮合刚度。问题在于一对齿从进入啮合到退出啮合的过程中啮合点沿齿廓不断移动轮齿的有效悬臂长度不停变化再加上重合度大于1时双齿啮合区和单齿啮合区交替出现参与承载的齿对数也在切换。于是综合啮合刚度就成了一条周期性波动的曲线。这条曲线直接以参数激励的形式进入齿轮动力学方程成为振动响应的主要激励源。工程上常见的齿轮箱噪声超标、齿面点蚀加速、轴承过早失效追根溯源往往都能在刚度曲线的波动幅度和跳变位置上找到线索。1.2 石川算法的定位工程精度和计算成本的最佳折中计算啮合刚度的方法有好几条路线有限元法最精细能反映齿根过渡曲线、轮辐结构等细节但建模耗时修一个参数就要重新划分网格韦伯-巴纳谢克法把载荷沿齿廓分布处理精度也不错但要处理积分核函数代码复杂一些石川算法的思路是做最激进的简化——把轮齿当成一个变截面悬臂梁用材料力学公式把弯曲、剪切、压缩和齿基变形逐项算出来然后叠加。很多人担心石川算法精度不够。实际对比下来在标准直齿轮的正常工作范围内石川法算出的齿根弯曲应力偏大10%~20%但时变刚度的趋势和量级完全正确用于动力学响应和故障特征分析绰绰有余。而且它计算一个啮合周期的刚度只要毫秒级做参数扫描、优化设计、批量故障仿真非常舒服。下面这张表把几类方法的定位理了一下方法模型基础相对计算量适用场景有限元法真实齿形实体大应力精细校核、齿根裂纹扩展模拟韦伯法齿廓分布载荷中齿廓修形设计、接触应力分布研究石川算法变截面悬臂梁极小时变刚度快速计算、动力学迭代、优化扫描我在实际项目中通常先用石川法做全工况扫描找到风险点后再用有限元对关键工况精确复算。这套组合拳效率非常高。2. 石川算法的力学图景一根变截面悬臂梁2.1 轮齿的悬臂梁简化与变形分解石川算法最核心的假设是单个轮齿可以看作固定在齿根圆上的悬臂梁载荷作用于齿廓啮合点方向沿啮合线法向。齿面形状从齿根到齿顶逐渐变薄所以梁的截面高度即齿厚沿径向变化这是变截面的含义。在这个模型下轮齿在法向载荷下的总柔度力与变形之比取倒数由五部分串联而成弯曲变形、剪切变形、轴向压缩变形、齿根基础弹性变形、以及两齿面接触处的赫兹弹性压扁。串联的意思是这些变形发生在同一条传力路径上总柔度等于各部分柔度之和物理上可以类比为弹簧组逐个叠放总柔度相加。2.2 五部分柔度的数学表达以主动轮为例任意半径 r 处的截面齿厚记为 s(r)截面惯性矩 I(r)b·s(r)³/12其中 b 为齿宽。设载荷作用点半径为 r_c齿根约束位置半径为 r_f工作压力角为 α_w沿径向的柔度积分可以写成弯曲柔度近似为 cos²α_w 乘上从 r_f 到 r_c 的 (r_c−r)²/(E·I(r)) 对 r 的积分这部分占据了轮齿弹性变形的绝大部分物理上就是悬臂梁弯矩产生的转角与挠度。剪切柔度1.2·(1ν)·cos²α_w 乘上对 1/(b·s(r)) 的积分并除以弹性模量 E系数1.2是矩形截面剪应力分布不均匀修正。压缩柔度sin²α_w 乘上对 1/s(r) 的积分再除以 (E·b)处理的是法向载荷沿径向分量把轮齿压短的效应。齿基柔度用一个工程经验公式估算比较常见的写法是 cos²α_w/(E·b) 乘上 ((l_f/u_f)² 2.4/(1ν) tan²α_w)其中 l_f 是载荷点到齿根的距离u_f 是齿根截面齿厚。赫兹接触柔度两齿面在啮合点发生线接触弹性压扁接触刚度取 k_hπ·E·b/(4(1−ν²))柔度直接取 1/k_h。从动轮同样要计算前四部分整套公式里把主动轮和从动轮的几何参数分别代进去最后五部分全部相加取倒数就得到一对齿在某个啮合位置的啮合刚度。2.3 为什么这样拆分是合理的我在初学的时候觉得这套公式很粗糙但后来发现它的每一步都有清晰物理含义。弯曲和剪切处理的是轮齿作为悬臂梁的主体变形压缩处理径向力分量齿基柔度弥补了刚性固定端假设带来的误差赫兹项则处理接触局部的非线性弹性。五部分相互独立又相互补充没有一个冗余项。更重要的是这种拆分方式让故障仿真变得很直观——想模拟齿根裂纹就削减裂纹截面的惯性矩想模拟点蚀剥落就降低赫兹接触刚度。这种哪里坏了改哪里的便利性是有限元模型很难比拟的。3. 开始编码前的几何准备让MATLAB认识这对齿轮3.1 齿轮基本参数与几何量计算先定义一对标准安装的直齿圆柱齿轮副。为了让结果有代表性我选了模数2.5、主动轮24齿、从动轮48齿、压力角20°、零变位的组合齿宽20mm材料弹性模量206000MPa泊松比0.3主动轮输入转矩100N·m。这套参数算出来的重合度在1.6~1.7之间单双齿啮合区间划分很典型。m 2.5; z1 24; z2 48; alpha_n 20/180*pi; x1 0; x2 0; b 20; E 206e3; nu 0.3; T1 100e3; % N*mm d1 m*z1; d2 m*z2; % 分度圆直径 db1 d1*cos(alpha_n); db2 d2*cos(alpha_n); % 基圆直径 da1 d1 2*(1x1)*m; da2 d2 2*(1x2)*m; % 齿顶圆直径 df1 d1 - 2*(1.25-x1)*m; df2 d2 - 2*(1.25-x2)*m; % 齿根圆直径 a0 (d1 d2)/2; % 标准中心距 alpha_w alpha_n; % 标准安装时工作压力角等于分度圆压力角这里有个细节容易踩坑如果齿轮副不是标准安装也就是实际中心距不等于标准中心距那么工作压力角必须由中心距反算公式是 acos(a0·cosα_n/a_w)其中 a_w 是实际中心距。很多工程案例里为了调整侧隙会取大一点的中心距这时候直接用 α_n 算出来的刚度误差比较明显。3.2 渐开线齿廓与任意半径齿厚石川算法需要知道齿廓上任意半径位置的齿厚。渐开线齿轮的齿厚公式并不复杂% 渐开线函数 inv (u) tan(u) - u; % 分度圆齿厚 s0 pi*m/2 2*x1*m*tan(alpha_n); % 任意半径 r 处的齿厚r rb/2 alpha_r acos(db1./(2*r)); s_r r .* (s0/(d1/2) 2*(inv(alpha_n) - inv(alpha_r)));理解这条公式的关键在于从分度圆到任意半径齿厚按渐开线展角差修正。半径越大齿厚越小半径小于基圆时齿面不再是渐开线公式里的反余弦会超出定义域。处理办法是让积分起点取 max(齿根圆半径, 基圆半径)把基圆以下的齿根过渡区当作刚性附着区。这会让计算结果偏刚一点但趋势完全正确后面我会专门说这个近似的影响。3.3 啮合线长度、基圆齿距与单双齿区间啮合几何是整个程序的地基。需要先算出实际啮合线长度、基圆齿距和重合度rho_a1 sqrt((da1/2)^2 - (db1/2)^2); % 主动轮齿顶处曲率半径 rho_a2 sqrt((da2/2)^2 - (db2/2)^2); % 从动轮齿顶处曲率半径 g_a rho_a1 rho_a2 - a0*sin(alpha_w); % 实际啮合线长度 p_b pi*m*cos(alpha_n); % 基圆齿距 epsilon g_a / p_b; % 端面重合度实际啮合线上的起始点是主动轮齿根与从动轮齿顶接触的位置终点是主动轮齿顶与从动轮齿根接触的位置。从起点沿啮合线走一个基圆齿距 p_b第二对齿才进入啮合。所以对一副重合度在1到2之间的直齿轮来说前半段和后半段各有一段双齿啮合区中间是单齿啮合区。双齿区长度 g_a−p_b单齿区长度 2p_b−g_a。判断当前啮合位置上参与啮合的齿对数最稳妥的写法是s_grid linspace(0, g_a, 200); for i 1:length(s_grid) s s_grid(i); if s p_b || s g_a - p_b n_contact 2; else n_contact 1; end end这段逻辑虽然简单但它是后面刚度叠加的基础。尤其要注意一个容易弄反的地方s0 对应的是这对齿刚开始啮合而不是节点位置。很多人习惯以节点为坐标原点导致后面 rho1、rho2 的符号搞错。我建议统一用从啮合起点出发的啮合线位移作为自变量整个程序从头到尾不要换坐标系这样最不容易乱。4. 核心代码实现从柔度积分到刚度曲线4.1 单齿对刚度函数石川算法里工作量最大的部分就是单齿对刚度计算函数。我把它封装成single_pair_stiffness输入是该齿对从进入啮合起的局部位移 s输出是当前啮合位置的单齿对啮合刚度。function k single_pair_stiffness(s, geom1, geom2, mat) % 1. 由啮合位移求两轮接触点曲率半径 rho_a2 sqrt((geom2.da/2)^2 - (geom2.db/2)^2); rho1_start geom1.g_a - rho_a2; % 啮合起点处主动轮曲率半径 rho2_start rho_a2; % 啮合起点处从动轮曲率半径齿顶处 rho1 rho1_start s; rho2 rho2_start - s; % 2. 接触点半径 r1c sqrt((geom1.db/2)^2 rho1^2); r2c sqrt((geom2.db/2)^2 rho2^2); % 3. 对两个轮子分别做梁积分 d1 bending_shear_axial(geom1, r1c, mat); d2 bending_shear_axial(geom2, r2c, mat); % 4. 齿基柔度 d_f1 foundation_flex(geom1, r1c, mat); d_f2 foundation_flex(geom2, r2c, mat); % 5. 赫兹接触柔度 kh pi*mat.E*geom.b / (4*(1-mat.nu^2)); d_h 1/kh; % 6. 总柔度取倒数得到刚度 k 1 / (d1 d2 d_f1 d_f2 d_h); end这里的关键点是 rho1_start 和 rho2_start 的方向。从啮合起点出发主动轮接触点沿啮合线向齿顶方向移动所以曲率半径增加从动轮接触点从齿顶向齿根移动曲率半径减小。我把所有符号统一到一个坐标系里避免了在两个齿轮之间反复切换方向带来的混乱。4.2 弯曲、剪切、轴向柔度的数值积分bending_shear_axial函数负责从齿根约束位置到接触点位置做离散积分。我用40个截面点做梯形积分收敛性完全足够function d_total bending_shear_axial(geom, r_contact, mat) r_root max(geom.df/2, geom.db/2); % 从基圆或齿根圆开始 rv linspace(r_root, r_contact, 40); alpha_r acos(geom.db./(2*rv)); s_r rv .* (geom.s0/(geom.d/2) 2*(tan(pi/6) - pi/6 ... - (tan(alpha_r) - alpha_r))); % 写成调用内部inv函数更清晰 I_r mat.b .* s_r.^3 / 12; A_r mat.b .* s_r; alpha_w geom.alpha_w; d_bend cos(alpha_w)^2 * trapz(rv, (r_contact - rv).^2 ./ (mat.E*I_r)); d_shear 1.2 * (1mat.nu) / mat.E * cos(alpha_w)^2 * trapz(rv, 1./A_r); d_axial sin(alpha_w)^2 / (mat.E*mat.b) * trapz(rv, 1./s_r); d_total d_bend d_shear d_axial; end从结果占比来看弯曲柔度大约占60%~70%剪切次之轴向压缩通常不到5%。这也是为什么很多简化算法只算弯曲和剪切也能凑合用但要做到和实测趋势一致剪切和齿基项还是不能省。40个截面点听起来很少但梯形积分对这个光滑函数收敛已经足够我之前试过200个点刚度变化不到0.5%白浪费计算时间。齿基柔度的经验公式直接按下面代码写function d_f foundation_flex(geom, r_contact, mat) r_root max(geom.df/2, geom.db/2); l_f r_contact - r_root; u_f geom.df; % 用齿根圆处近似的等效齿宽截面 alpha_w geom.alpha_w; d_f cos(alpha_w)^2/(mat.E*geom.b) * ... ( (l_f/u_f)^2 2.4/(1mat.nu) tan(alpha_w)^2 ); end4.3 啮合周期循环与双齿区叠加单齿对刚度算好以后时变啮合刚度就是把同一时刻参与啮合的所有齿对刚度叠加起来。对重合度小于2的直齿轮规则清晰第一对齿从 s0 到 sg_a 全程参与第二对齿从 sp_b 开始参与到 sg_a 退出。于是总刚度按下面的循环计算k_total zeros(size(s_grid)); for i 1:length(s_grid) s s_grid(i); k1 single_pair_stiffness(s, geom1, geom2, mat); % 第一对齿 if s p_b k2 single_pair_stiffness(s - p_b, geom1, geom2, mat); % 第二对齿 k_total(i) k1 k2; else k_total(i) k1; end end注意第二对齿传入的是 s−p_b也就是它从自己进入啮合点开始的相对位移。这是最容易错的地方我第一版程序直接用 s 代入第二对齿结果双齿区刚度严重偏大检查了很久才发现是相位没对齐。叠加上去以后你会看到典型的时变刚度曲线两端双齿区刚度高中间单齿区刚度低整体呈现两边高中间低的形状。如果重合度再大一点比如1.8中间单齿区很短曲线更像一个浅V形。4.4 一个值得注意的代码边界问题前面提到积分起点取 max(齿根圆半径, 基圆半径)这个选择背后有代价。如果齿轮齿数较少比如z117以下齿根圆可能明显小于基圆石川法把基圆到齿根之间的过渡段当成刚性体相当于人为缩短了悬臂梁长度。结果就是齿数少时算出的刚度偏大。我在做少齿数齿轮时会在齿根过渡区补一段直线齿廓近似把齿厚从基圆处线性收窄到齿根圆处这样算出的刚度更接近有限元结果。对常规齿数z20来说这个修正量不到5%不做也没关系。5. 齿面接触变形与载荷分担刚度之外的第二张图5.1 赫兹接触变形的物理意义很多人算完刚度就停了但齿面接触变形量同样重要。所谓接触变形就是两齿面在啮合点处因赫兹接触产生的局部弹性压扁量。它虽然只占齿轮总变形的一小部分却直接关系到齿面胶合、点蚀的判别也影响实际载荷在齿对之间的分配。赫兹接触柔度在石川框架里是单独一项数值大小取决于材料参数和齿宽与啮合位置基本无关。按前面公式计算本案例的接触刚度 k_h 大约在3.5×10⁶ N/mm量级比单齿对综合刚度大约10⁵ N/mm量级大一个数量级以上。所以接触柔度在总柔度里占比大约10%~15%不是主导项但在载荷分担和故障模拟里不能省略。5.2 载荷分担率双齿区里力是怎么分的啮合过程中法向总载荷在参与啮合的齿对之间分配分配比例近似按各齿对刚度比例。单齿区当然全部载荷由一对齿承担双齿区则由两对齿按刚度分担。写成代码很直接F_total T1/(d1/2); for i 1:length(s_grid) if s_grid(i) p_b || s_grid(i) g_a - p_b s2_local s_grid(i) - p_b; k1 single_pair_stiffness(s_grid(i), geom1, geom2, mat); k2 single_pair_stiffness(s2_local, geom1, geom2, mat); F1 F_total * k1 / (k1 k2); F2 F_total - F1; delta_h1 F1 / kh; delta_h2 F2 / kh; else F1 F_total; delta_h1 F1 / kh; end end画出接触变形曲线后你会发现单齿区的接触变形明显大于双齿区因为全部载荷压在一对齿上而双齿区两对齿分摊载荷虽然各自刚度接近但每对齿的接触变形显著下降。这条曲线和刚度曲线的形状正好相反是解释齿面载荷波动最直观的材料。更精细的动力学分析里载荷分担率还会受到误差激励和振动惯性的影响静态刚度分配只是一个基准但绝大多数故障诊断和强度校核场景下这个基准已经足够定位问题。5.3 与ISO标准算法的对比经验ISO 6336标准里用啮合刚度评价齿轮承载能力单位是N/(mm·μm)也就是每单位齿宽每微米变形对应的载荷。我做过一组对比用石川法算出总变形后换算成单位齿宽刚度和ISO算法比石川法结果通常低10%~15%主要是齿根过渡区简化和齿基柔度经验公式带来的偏差。这个偏差在工程判断里完全可以接受但要注意一致性——如果你用石川法做优化对比所有方案都用同一套算法相对趋势是可靠的如果要和标准许用值直接对比建议做一次系数校准。6. 看懂振动响应图表从刚度和变形到动力学特征6.1 时变刚度曲线和接触变形曲线怎么读刚度和变形曲线不能只看形状要会读里面的工程信息。时变刚度曲线最典型的特征就是两端高、中间低的凹陷。这个凹陷深度也就是单齿区刚度与双齿区刚度的差值直接决定了齿轮啮合冲击的能量大小。差值越大啮入啮出时刻的冲击越剧烈振动幅值越高。设计时如果发现单齿区刚度过低通常要增加齿宽或增大模数而不是盲目修改变位系数。接触变形曲线的波峰值对应的是最不利承载位置。如果这个位置恰好落在齿顶或者齿根说明齿廓修形参数可能需要调整。我在诊断那个减速机异响案例时就是看到接触变形曲线的峰值位置偏移到齿顶附近才判断齿顶修缘量不足和现场齿面斑点互相印证。6.2 用Newmark-β法搭一个单自由度齿轮动力学模型有了时变刚度曲线可以进一步做振动响应。最简单的模型是单自由度齿轮副扭转振动模型沿啮合线方向的相对位移作为广义坐标动力学方程为m_e·x c·x k(t)·x F_m e(t)·k(t)其中 m_e 是等效质量对直齿轮可以由两轮转动惯量折算到啮合线方向c 是啮合阻尼通常取刚度临界阻尼的3%~10%k(t) 就是前面算出的时变啮合刚度F_m 是平均法向载荷e(t) 是综合啮合误差可以取定值也可以按谐波模拟。用Newmark-β法做时间步进MATLAB里核心循环不复杂dt 2e-6; t_end 0.2; tt 0:dt:t_end; x0 0; v0 0; fm F_total * cos(alpha_w); for n 1:length(tt)-1 t_n tt(n); s_phase mod(t_n * f_m, 1); % 归一化啮合相位 k_n interp1(s_grid/g_a, k_total, s_phase, linear, extrap); [x1, v1, a1] newmark_step(k_n, c, m_eq, fm, x0, v0, dt, beta, gamma); x0 x1; v0 v1; end这里最关键的一步是interp1把空间域的刚度曲线转换成时间域的激励。啮合频率 f_m z1·f_rot其中 f_rot 是主动轮轴频。转速一变刚度激励的基频就变系统的响应也随之改变。做升速试验仿真时只需在时间循环里按斜率更新 f_rot 即可程序不用大改。6.3 故障齿轮的刚度特征与频谱边频石川算法的几何参数化特性让故障仿真非常简单。齿根裂纹是最常见的故障模式在计算弯曲柔度时把裂纹位置所在截面的齿厚从 s_r 减小到 s_r − 2a其中 a 是裂纹深度沿齿宽方向的投影长度。截面惯性矩随之下降弯曲柔度增加体现在刚度曲线上就是该齿对应啮合相位出现一个局部凹陷。裂纹越深凹陷越明显退出啮合时刚度恢复的斜率越陡。齿面点蚀对刚度的影响路径不同。点蚀改变的是接触表面质量在模型里可以近似为赫兹接触刚度的局部折减。把点蚀区域对应的啮合相位范围内的 k_h 乘以一个0.6~0.9的衰减系数就能模拟出刚度曲线上的局部凹坑。故障齿轮做频谱分析时最值得关注的不是啮合频率本身的幅值而是它两侧的边频带。健康齿轮的频谱在啮合频率 f_m 处有离散谱线伴随少量均匀分布的齿轮固有频率响应裂纹齿轮由于每转一圈产生一次局部刚度损失会在 f_m 两侧以轴频 f_rot 为间隔出现一系列边频幅值随故障程度增加。如果把时域信号做阶次谱边频位置对应 z1±1、z1±2阶这个特征在诊断现场非常好用。我用这个流程处理过含裂纹齿轮的数据仿真频谱上的边频带分布和实测振动信号高度一致差别主要在幅值比例上——因为实际齿轮还有轴系不对中、轴承故障等其他激励叠加。这反过来说明把石川刚度作为激励源放进动力学模型做故障机理研究路径是通顺的。7. 实测中的几个坑以及一点个人建议7.1 单位体系混乱是振动分析报错的第一大来源这个坑我必须放在第一个说。齿轮计算中长度用毫米、压力用兆帕、力的单位用牛顿是最顺手的一套组合但柔度积分公式里的 E·I 计算、赫兹接触刚度的单位换算稍有混用就会出现量级离谱的结果。我见过有人把 E 用成 2.06×10¹¹ Pa其他量却用毫米算出来的柔度直接小了1000倍刚度曲线画出来跟心电图似的。建议在代码开头统一用mm、N、MPa、N/mm这套单位并在注释里写明每个输入变量的单位和量级。7.2 重合度接近整数时刚度曲线为什么会抖当齿轮副参数使重合度接近2.0时单齿区几乎消失双齿区覆盖整个啮合周期刚度曲线不再有明显的大幅凹陷而是在平均刚度附近叠加小幅波动。这时候如果离散点数太少比如只用50个点曲线会出现锯齿状抖动看起来像数值不稳定实际上只是采样分辨率不够。把 s_grid 的点数提高到500以上或者改用非均匀采样就能得到光滑的波形。另一个容易误判的情况是重合度刚好跨越整数时齿对数的切换位置会产生刚度突变跳点这个不是bug而是真实的冲击激励分析时反而要重点标注。7.3 建议先用文献曲线验证再谈精度最后说点实在的。石川算法的代码写起来不难难的是确认你的实现是对的。我的建议是不要一上来就跑自己的工程参数先用一组教科书上的标准齿轮参数比如模数4、z120、z230、压力角20°、零变位把算出的最小单齿刚度、最大双齿刚度和文献里的参考曲线对比。量级一致、趋势吻合再切换到自己的数据。这一步能省掉后面排查程序bug的大量时间。还有个小技巧做故障对比研究时固定一套基准参数只改变裂纹深度或点蚀位置观察刚度和频谱响应的相对变化而不是纠结绝对精度。相对变化趋势才是故障特征提取最有价值的信息。我这些年用石川算法踩过的坑基本都集中在几何边界处理和单位换算上物理模型本身反而很少出问题。这套方法虽然不如有限元精致但在工程快速评估和故障机理研究里它的小快灵优势是实打实的。拿到一组齿轮参数几分钟内就能从刚度曲线看到啮合冲击强弱从接触变形曲线判断齿面载荷分布再从动力学频谱定位激励特征整个链条走下来对一个齿轮系统的状态心里基本就有数了。