简介这份MATLAB脚本面向齿轮动力学与机械振动方向的研究生、工程师及科研人员用于求解直齿轮啮合过程中随时间动态变化的啮合刚度。核心算法基于势能法通过建立轮齿变形能与外载荷之间的关系结合齿面接触几何、材料弹性特性及啮合位置变化逐点计算不同时刻的刚度数值可生成单对齿啮入至啮出全过程的刚度曲线输出时间序列对应的刚度数组直接服务于齿轮系统动力学建模、振动响应分析与故障特征提取。资源包共6个文件约42KB包含m主计算脚本、py辅助脚本、txt说明文档、inscode工程配置及结果示意图等结构清晰输入参数涵盖模数、齿数、压力角、齿宽、弹性模量、泊松比与转速等基础齿轮参数无需额外工具箱在MATLAB基础平台即可运行。目前已有36人学习下载适合需要快速获取时变啮合刚度数据、开展后续仿真与验证工作的读者参考使用。1. 直齿轮啮合刚度随时间变化为什么势能法值得你花一个下午跑通如果你正在做齿轮箱动力学仿真大概率遇到过这个尴尬振动响应算出来和实测对不上频谱里该出现的啮合频率边带要么位置偏了要么幅值差一大截。排查一圈轴承、箱体、联轴器都没问题最后发现根子在激励源——啮合刚度被当成了常数。真实齿轮在啮合过程中参与承载的齿对数在单齿对和双齿对之间来回切换刚度是随时间周期性波动的这个波动才是振动里那些边带和调制的来源。势能法Potential Energy Method就是算这个时变刚度最常用的一条路。它的思路很直接把一对轮齿看成一端固定的变截面悬臂梁分别算出弯曲势能、剪切势能、轴向压缩势能和接触势能再对各自求偏导得到对应的柔度分量最后取倒数并叠加就得到单齿对的啮合刚度。相比有限元接触分析势能法不需要建模、不需要求解器、不需要几小时的网格划分一个 MATLAB 脚本几十行就能跑出整条刚度曲线参数改起来也快。适合做动力学建模需要快速拿到刚度激励的从业者也适合想搞明白刚度曲线每个拐点从哪来的新手。下面把我自己跑这套东西的完整路径拆开讲。2. 势能法算啮合刚度的力学模型与公式落地2.1 把轮齿拆成悬臂梁四个柔度分量从哪来势能法的核心假设是轮齿的变形由四部分组成每一部分对应一种能量存储方式。弯曲势能来自齿在法向力作用下像悬臂梁一样弯剪切势能来自截面剪切变形轴向压缩势能来自压力沿齿高方向的压缩接触势能来自两齿面接触区的局部变形。前三个是轮齿本体的第四个是接触的。对渐开线直齿轮齿廓上任意位置的法向力可以分解为沿齿高和垂直齿高两个方向。设啮合点作用力为 F作用线与齿中心线夹角为 α那么弯曲、剪切、压缩三个柔度分量可以写成沿齿高积分的解析式。接触柔度用赫兹接触理论只和材料弹性模量、泊松比、齿宽有关与啮合位置无关。这里有个容易翻车的地方积分上下限。齿根圆到啮合点之间的齿廓才是参与变形的部分积分要从齿根圆半径积到当前啮合点半径。很多人直接把整个齿高积分结果刚度偏小一大截。我一般会先把齿廓离散成若干点用数值积分代替解析积分这样齿形修形、齿根过渡曲线都能直接塞进去通用性更好。2.2 单齿对刚度与重合度叠加从一条曲线到整条时变曲线单齿对啮合刚度 k_single 是四个柔度分量之和的倒数k_single 1 / (δ_b δ_s δ_a δ_h)其中 δ_b、δ_s、δ_a 分别是弯曲、剪切、压缩柔度δ_h 是接触柔度。注意这里每个柔度都是啮合点位置的函数所以 k_single 随啮合点移动而变化通常在齿根附近刚度大、齿顶附近刚度小。但真实齿轮不是只有一对齿在啮合。重合度 ε 大于 1 时会有一段区间两对齿同时承载。总刚度不是简单相加而是两对齿的柔度先并联再取倒数k_total 1 / (1/k_pair1 1/k_pair2)这里 k_pair1 和 k_pair2 分别是两对齿各自的单齿对刚度但要注意它们对应的啮合位置不同一个在齿根附近一个在齿顶附近。所以整条时变刚度曲线的形状本质上是单齿对刚度曲线按重合度平移叠加的结果。重合度越大双齿啮合区越长刚度波动越平缓。这也是为什么斜齿轮比直齿轮振动小——斜齿轮重合度大刚度波动被抹平了。2.3 参数表先对齐这 8 个量再动手在写代码之前把下面这张表填好。参数不对后面全白算。参数符号典型取值说明模数m2 mm决定齿高和齿厚齿数z1, z220, 40决定传动比和重合度压力角α20°标准值影响齿廓形状齿宽B20 mm直接影响接触柔度弹性模量E206 GPa钢齿轮常用值泊松比ν0.3钢齿轮常用值转速n1500 rpm决定时间轴映射重合度ε1.6 左右由齿数压力角算出重合度可以用经验公式先估ε ≈ [z1(tan α_a1 - tan α) z2(tan α_a2 - tan α)] / (2π)其中 α_a 是齿顶压力角。这个值决定了双齿啮合区占整个啮合周期的比例是后面叠加逻辑的关键输入。3. MATLAB 脚本实现从齿廓离散到刚度曲线输出3.1 齿廓离散与积分变量准备先把渐开线齿廓离散成点同时算出每个点对应的半径、齿厚、压力角。这一步是整个脚本的地基。% 基本参数 m 2e-3; % 模数 (m) z1 20; z2 40; % 齿数 alpha 20*pi/180; % 压力角 (rad) B 20e-3; % 齿宽 (m) E 206e9; % 弹性模量 (Pa) nu 0.3; % 泊松比 % 分度圆、基圆、齿顶圆、齿根圆半径 r1 m*z1/2; r2 m*z2/2; rb1 r1*cos(alpha); rb2 r2*cos(alpha); ra1 r1 m; ra2 r2 m; rf1 r1 - 1.25*m; rf2 r2 - 1.25*m; % 齿廓离散从齿根圆到齿顶圆 N 200; r_profile linspace(rf1, ra1, N); % 渐开线极角 theta sqrt((r_profile/rb1).^2 - 1) - acos(rb1./r_profile); % 齿厚半角分度圆处齿厚一半对应的角度 s_half pi/(2*z1); phi s_half theta; % 齿廓上每点的半齿厚角 y r_profile .* sin(phi); % 齿厚方向坐标这段代码的逻辑是渐开线上任意一点的极角由基圆半径和该点半径决定齿厚半角等于分度圆半齿厚角加上渐开线展开角。y 就是该点处齿厚的一半后面算截面惯性矩和截面积要用。参数 N 控制离散精度200 个点对一般模数够用模数特别小或者修形复杂时可以加到 500。3.2 四个柔度分量的数值积分有了齿廓离散点就可以对每个啮合位置算四个柔度。下面这段是核心计算。% 啮合点位置沿啮合线从齿根到齿顶 x_mesh linspace(rf1, ra1, N); delta_b zeros(1,N); delta_s zeros(1,N); delta_a zeros(1,N); for i 1:N % 当前啮合点半径 rx x_mesh(i); % 当前点对应的齿厚半角和力臂 theta_x sqrt((rx/rb1)^2 - 1) - acos(rb1/rx); phi_x s_half theta_x; yx rx * sin(phi_x); % 力臂啮合点到齿根的距离在力方向上的投影 L rx * cos(phi_x) - rf1 * cos(s_half); % 对齿根到当前点之间的齿廓积分 idx r_profile rx; r_int r_profile(idx); y_int y(idx); % 截面惯性矩和面积 I_int (2*y_int).^3 * B / 12; A_int 2*y_int * B; % 弯曲柔度积分 delta_b(i) trapz(r_int, L.^2 ./ (E * I_int)); % 剪切柔度积分引入剪切修正系数 kappa 1.2; delta_s(i) trapz(r_int, kappa ./ (E * A_int)); % 压缩柔度积分 delta_a(i) trapz(r_int, 1 ./ (E * A_int)); end % 接触柔度赫兹接触与位置无关 delta_h 4*(1-nu^2) / (pi * E * B);逻辑说明弯曲柔度积分里 L 是力臂随啮合点变化I_int 是截面惯性矩随齿厚变化。剪切柔度里 kappa 是剪切修正系数矩形截面取 1.2。压缩柔度最简单只和截面积有关。接触柔度用赫兹公式对钢-钢接触这个值通常在 1e-8 量级比前三个小一到两个数量级但齿顶附近弯曲柔度变大时它就不能忽略了。参数说明trapz 是梯形数值积分比 sum 精度高。如果齿廓点不够密积分误差会累积表现为刚度曲线毛刺。我一般会检查刚度曲线是否光滑有毛刺就加 N。3.3 单齿对刚度与双齿叠加逻辑四个柔度算完单齿对刚度就是倒数。然后按重合度做平移叠加。% 单齿对刚度 k_single 1 ./ (delta_b delta_s delta_a delta_h); % 重合度估算 alpha_a1 acos(rb1/ra1); alpha_a2 acos(rb2/ra2); eps_a (z1*(tan(alpha_a1)-tan(alpha)) ... z2*(tan(alpha_a2)-tan(alpha))) / (2*pi); % 啮合周期一个齿距对应的时间 n 1500; % 转速 rpm omega 2*pi*n/60; % 角速度 rad/s T_mesh 2*pi/(omega*z1); % 啮合周期 (s) % 时间轴 t linspace(0, 2*T_mesh, 1000); % 双齿叠加第二对齿的刚度曲线是第一对的平移 shift round(eps_a * N); % 平移点数 k_total zeros(size(t)); for j 1:length(t) % 映射到单齿对曲线上的位置 pos mod(t(j)/T_mesh, 1) * N; idx1 max(1, min(N, round(pos)1)); k1 k_single(idx1); % 第二对齿位置偏移一个重合度 pos2 pos eps_a*N; if pos2 N idx2 round(pos2)1; k2 k_single(max(1,min(N,idx2))); k_total(j) 1/(1/k1 1/k2); else k_total(j) k1; end end plot(t*1e3, k_total/1e6, LineWidth, 1.5); xlabel(时间 (ms)); ylabel(啮合刚度 (MN/m)); grid on;逻辑说明单齿对刚度曲线是一条随啮合位置变化的曲线。双齿啮合时第二对齿的啮合位置比第一对超前一个重合度对应的弧长所以把单齿对曲线平移 eps_a*N 个点就是第二对齿的刚度。两对齿并联柔度相加再取倒数。当第二对齿还没进入啮合或已经退出时只有单齿对承载直接取 k1。参数说明shift 是平移点数由重合度决定。eps_a 算出来通常在 1.5 到 1.8 之间对应双齿啮合区占 50% 到 80%。时间轴取两个啮合周期是为了看周期性。刚度单位用 MN/m 是因为钢齿轮单齿对刚度典型值在 1e8 到 1e9 N/m 量级。4. 避坑与排查刚度曲线不对时先查这 5 处4.1 刚度值整体偏大或偏小一个数量级现象跑出来的刚度峰值只有 1e6 N/m 或者高达 1e11 N/m和文献里的 1e8 到 1e9 对不上。原因单位不统一。模数用了毫米但弹性模量用了帕或者齿宽用了毫米但没转成米。势能法里长度单位必须统一到米力单位统一到牛顿。解决在脚本开头把所有几何量统一转成米E 用帕。跑完先看单齿对刚度峰值钢齿轮齿宽 20 mm 时应该在 5e8 N/m 左右。差一个数量级就是单位问题差两个数量级检查积分上下限。4.2 刚度曲线出现周期性毛刺现象曲线整体趋势对但上面叠加了高频锯齿。原因齿廓离散点太少数值积分在齿根附近截面变化剧烈的地方误差大。或者平移叠加时索引取整引入了跳变。解决把 N 从 200 加到 500毛刺通常会消失。如果还有检查平移逻辑里 idx2 的计算用线性插值代替最近点取整。4.3 双齿啮合区刚度反而比单齿区低现象叠加后双齿区的刚度比单齿区还小违反直觉。原因并联公式写反了。两对齿并联是柔度相加再取倒数不是刚度直接相加。如果写成 k1k2双齿区会异常大如果写成 1/(k1k2)双齿区会异常小。解决确认公式是 1/(1/k1 1/k2)。这个式子等价于 k1*k2/(k1k2)永远小于 min(k1,k2)但大于单对齿柔度倒数的一半。4.4 重合度算出来小于 1现象eps_a 算出 0.8 之类导致没有双齿啮合区。原因齿顶圆半径算错或者压力角用错。标准直齿轮重合度一般大于 1.2小于 1 说明齿顶高系数或者压力角输入有误。解决检查 ra r m 是否用了正确的齿顶高系数标准是 1。压力角确认是 20 度不是 25 度。如果确实设计成小重合度那单齿啮合区就是全部叠加逻辑要相应简化。4.5 时间轴和啮合周期对不上现象刚度曲线周期和理论啮合周期不一致或者一个周期里出现了两次波动。原因啮合周期公式用错。啮合周期是转过一个齿距的时间等于 2π/(ω*z)不是 2π/ω。z 用主动轮齿数。解决确认 T_mesh 2pi/(omegaz1)。如果曲线一个周期里有两个波峰说明时间轴取长了或者叠加逻辑里平移量算错导致重复叠加。5. 进阶把刚度曲线喂给动力学模型并验证5.1 从刚度曲线到动力学激励拿到 k_total(t) 只是第一步真正要用它得把它变成动力学方程里的时变系数。常见做法是把刚度曲线做傅里叶展开取前几阶谐波作为激励% 对刚度曲线做 FFT K fft(k_total - mean(k_total)); K_mag abs(K)/length(t)*2; f (0:length(t)-1)/(t(end)-t(1)); % 取前 5 阶谐波 harmonics 1:5; K_harm K_mag(harmonics); fprintf(前5阶刚度谐波幅值 (MN/m):\n); disp(K_harm/1e6);逻辑说明减去均值后做 FFT得到各阶谐波幅值。啮合频率处的幅值就是刚度波动的主要激励。前 5 阶通常够用更高阶幅值衰减很快。这些谐波幅值直接对应动力学方程里刚度矩阵的时变部分。参数说明采样点数要够至少覆盖 10 个啮合周期。如果只取 2 个周期做 FFT频率分辨率太差谐波幅值不准。我一般会跑 20 个周期再截取稳态段做分析。5.2 用能量守恒做一次自检势能法算出来的刚度可以用能量法反推变形量做自检。在啮合点施加一个已知法向力 F用柔度算出总变形 δ F/k然后算变形能 U 0.5Fδ。这个能量应该等于四个势能分量之和。如果对不上说明柔度计算有遗漏。% 自检取齿顶附近一个啮合点 F_test 1000; % N idx_test round(0.8*N); delta_total delta_b(idx_test) delta_s(idx_test) ... delta_a(idx_test) delta_h; k_test 1/delta_total; U_deform 0.5 * F_test^2 * delta_total; % 四个分量能量之和 U_sum 0.5*F_test^2 * (delta_b(idx_test) delta_s(idx_test) ... delta_a(idx_test) delta_h); fprintf(变形能: %.6f J, 分量和: %.6f J\n, U_deform, U_sum);这个自检看起来是恒等式但实际写代码时如果某个柔度分量漏乘了系数或者积分上下限写错这里就会暴露。我习惯每改一次参数就跑一次这个检查比盯着曲线看靠谱。5.3 一个容易被忽略的细节齿根过渡曲线标准渐开线公式只描述工作齿廓齿根到渐开线起点之间有一段过渡曲线。这段曲线对弯曲柔度贡献不小尤其是齿根附近。如果直接用渐开线从齿根圆开始积分会高估齿根刚度。常见做法是用齿条刀具圆角包络出过渡曲线或者用一段圆弧近似。我一般用圆弧近似半径取 0.38m对刚度峰值影响在 5% 以内但曲线形状会更接近实测。跑完这套脚本你会得到一条随啮合位置变化的刚度曲线能直接喂给齿轮动力学模型。参数改起来快试不同齿数、模数、齿宽对刚度波动的影响一个下午能跑几十组。真正花时间的不是写代码是搞清楚每个参数背后的物理意义以及曲线不对时知道去哪找原因。希望帮到你。本文还有配套的精品资源点击获取 
长河君 · 2026/10/18 16:20:57 · 阅读 1.2k 势能法计算直齿轮时变啮合刚度：MATLAB实现与避坑指南
简介这份MATLAB脚本面向齿轮动力学与机械振动方向的研究生、工程师及科研人员用于求解直齿轮啮合过程中随时间动态变化的啮合刚度。核心算法基于势能法通过建立轮齿变形能与外载荷之间的关系结合齿面接触几何、材料弹性特性及啮合位置变化逐点计算不同时刻的刚度数值可生成单对齿啮入至啮出全过程的刚度曲线输出时间序列对应的刚度数组直接服务于齿轮系统动力学建模、振动响应分析与故障特征提取。资源包共6个文件约42KB包含m主计算脚本、py辅助脚本、txt说明文档、inscode工程配置及结果示意图等结构清晰输入参数涵盖模数、齿数、压力角、齿宽、弹性模量、泊松比与转速等基础齿轮参数无需额外工具箱在MATLAB基础平台即可运行。目前已有36人学习下载适合需要快速获取时变啮合刚度数据、开展后续仿真与验证工作的读者参考使用。1. 直齿轮啮合刚度随时间变化为什么势能法值得你花一个下午跑通如果你正在做齿轮箱动力学仿真大概率遇到过这个尴尬振动响应算出来和实测对不上频谱里该出现的啮合频率边带要么位置偏了要么幅值差一大截。排查一圈轴承、箱体、联轴器都没问题最后发现根子在激励源——啮合刚度被当成了常数。真实齿轮在啮合过程中参与承载的齿对数在单齿对和双齿对之间来回切换刚度是随时间周期性波动的这个波动才是振动里那些边带和调制的来源。势能法Potential Energy Method就是算这个时变刚度最常用的一条路。它的思路很直接把一对轮齿看成一端固定的变截面悬臂梁分别算出弯曲势能、剪切势能、轴向压缩势能和接触势能再对各自求偏导得到对应的柔度分量最后取倒数并叠加就得到单齿对的啮合刚度。相比有限元接触分析势能法不需要建模、不需要求解器、不需要几小时的网格划分一个 MATLAB 脚本几十行就能跑出整条刚度曲线参数改起来也快。适合做动力学建模需要快速拿到刚度激励的从业者也适合想搞明白刚度曲线每个拐点从哪来的新手。下面把我自己跑这套东西的完整路径拆开讲。2. 势能法算啮合刚度的力学模型与公式落地2.1 把轮齿拆成悬臂梁四个柔度分量从哪来势能法的核心假设是轮齿的变形由四部分组成每一部分对应一种能量存储方式。弯曲势能来自齿在法向力作用下像悬臂梁一样弯剪切势能来自截面剪切变形轴向压缩势能来自压力沿齿高方向的压缩接触势能来自两齿面接触区的局部变形。前三个是轮齿本体的第四个是接触的。对渐开线直齿轮齿廓上任意位置的法向力可以分解为沿齿高和垂直齿高两个方向。设啮合点作用力为 F作用线与齿中心线夹角为 α那么弯曲、剪切、压缩三个柔度分量可以写成沿齿高积分的解析式。接触柔度用赫兹接触理论只和材料弹性模量、泊松比、齿宽有关与啮合位置无关。这里有个容易翻车的地方积分上下限。齿根圆到啮合点之间的齿廓才是参与变形的部分积分要从齿根圆半径积到当前啮合点半径。很多人直接把整个齿高积分结果刚度偏小一大截。我一般会先把齿廓离散成若干点用数值积分代替解析积分这样齿形修形、齿根过渡曲线都能直接塞进去通用性更好。2.2 单齿对刚度与重合度叠加从一条曲线到整条时变曲线单齿对啮合刚度 k_single 是四个柔度分量之和的倒数k_single 1 / (δ_b δ_s δ_a δ_h)其中 δ_b、δ_s、δ_a 分别是弯曲、剪切、压缩柔度δ_h 是接触柔度。注意这里每个柔度都是啮合点位置的函数所以 k_single 随啮合点移动而变化通常在齿根附近刚度大、齿顶附近刚度小。但真实齿轮不是只有一对齿在啮合。重合度 ε 大于 1 时会有一段区间两对齿同时承载。总刚度不是简单相加而是两对齿的柔度先并联再取倒数k_total 1 / (1/k_pair1 1/k_pair2)这里 k_pair1 和 k_pair2 分别是两对齿各自的单齿对刚度但要注意它们对应的啮合位置不同一个在齿根附近一个在齿顶附近。所以整条时变刚度曲线的形状本质上是单齿对刚度曲线按重合度平移叠加的结果。重合度越大双齿啮合区越长刚度波动越平缓。这也是为什么斜齿轮比直齿轮振动小——斜齿轮重合度大刚度波动被抹平了。2.3 参数表先对齐这 8 个量再动手在写代码之前把下面这张表填好。参数不对后面全白算。参数符号典型取值说明模数m2 mm决定齿高和齿厚齿数z1, z220, 40决定传动比和重合度压力角α20°标准值影响齿廓形状齿宽B20 mm直接影响接触柔度弹性模量E206 GPa钢齿轮常用值泊松比ν0.3钢齿轮常用值转速n1500 rpm决定时间轴映射重合度ε1.6 左右由齿数压力角算出重合度可以用经验公式先估ε ≈ [z1(tan α_a1 - tan α) z2(tan α_a2 - tan α)] / (2π)其中 α_a 是齿顶压力角。这个值决定了双齿啮合区占整个啮合周期的比例是后面叠加逻辑的关键输入。3. MATLAB 脚本实现从齿廓离散到刚度曲线输出3.1 齿廓离散与积分变量准备先把渐开线齿廓离散成点同时算出每个点对应的半径、齿厚、压力角。这一步是整个脚本的地基。% 基本参数 m 2e-3; % 模数 (m) z1 20; z2 40; % 齿数 alpha 20*pi/180; % 压力角 (rad) B 20e-3; % 齿宽 (m) E 206e9; % 弹性模量 (Pa) nu 0.3; % 泊松比 % 分度圆、基圆、齿顶圆、齿根圆半径 r1 m*z1/2; r2 m*z2/2; rb1 r1*cos(alpha); rb2 r2*cos(alpha); ra1 r1 m; ra2 r2 m; rf1 r1 - 1.25*m; rf2 r2 - 1.25*m; % 齿廓离散从齿根圆到齿顶圆 N 200; r_profile linspace(rf1, ra1, N); % 渐开线极角 theta sqrt((r_profile/rb1).^2 - 1) - acos(rb1./r_profile); % 齿厚半角分度圆处齿厚一半对应的角度 s_half pi/(2*z1); phi s_half theta; % 齿廓上每点的半齿厚角 y r_profile .* sin(phi); % 齿厚方向坐标这段代码的逻辑是渐开线上任意一点的极角由基圆半径和该点半径决定齿厚半角等于分度圆半齿厚角加上渐开线展开角。y 就是该点处齿厚的一半后面算截面惯性矩和截面积要用。参数 N 控制离散精度200 个点对一般模数够用模数特别小或者修形复杂时可以加到 500。3.2 四个柔度分量的数值积分有了齿廓离散点就可以对每个啮合位置算四个柔度。下面这段是核心计算。% 啮合点位置沿啮合线从齿根到齿顶 x_mesh linspace(rf1, ra1, N); delta_b zeros(1,N); delta_s zeros(1,N); delta_a zeros(1,N); for i 1:N % 当前啮合点半径 rx x_mesh(i); % 当前点对应的齿厚半角和力臂 theta_x sqrt((rx/rb1)^2 - 1) - acos(rb1/rx); phi_x s_half theta_x; yx rx * sin(phi_x); % 力臂啮合点到齿根的距离在力方向上的投影 L rx * cos(phi_x) - rf1 * cos(s_half); % 对齿根到当前点之间的齿廓积分 idx r_profile rx; r_int r_profile(idx); y_int y(idx); % 截面惯性矩和面积 I_int (2*y_int).^3 * B / 12; A_int 2*y_int * B; % 弯曲柔度积分 delta_b(i) trapz(r_int, L.^2 ./ (E * I_int)); % 剪切柔度积分引入剪切修正系数 kappa 1.2; delta_s(i) trapz(r_int, kappa ./ (E * A_int)); % 压缩柔度积分 delta_a(i) trapz(r_int, 1 ./ (E * A_int)); end % 接触柔度赫兹接触与位置无关 delta_h 4*(1-nu^2) / (pi * E * B);逻辑说明弯曲柔度积分里 L 是力臂随啮合点变化I_int 是截面惯性矩随齿厚变化。剪切柔度里 kappa 是剪切修正系数矩形截面取 1.2。压缩柔度最简单只和截面积有关。接触柔度用赫兹公式对钢-钢接触这个值通常在 1e-8 量级比前三个小一到两个数量级但齿顶附近弯曲柔度变大时它就不能忽略了。参数说明trapz 是梯形数值积分比 sum 精度高。如果齿廓点不够密积分误差会累积表现为刚度曲线毛刺。我一般会检查刚度曲线是否光滑有毛刺就加 N。3.3 单齿对刚度与双齿叠加逻辑四个柔度算完单齿对刚度就是倒数。然后按重合度做平移叠加。% 单齿对刚度 k_single 1 ./ (delta_b delta_s delta_a delta_h); % 重合度估算 alpha_a1 acos(rb1/ra1); alpha_a2 acos(rb2/ra2); eps_a (z1*(tan(alpha_a1)-tan(alpha)) ... z2*(tan(alpha_a2)-tan(alpha))) / (2*pi); % 啮合周期一个齿距对应的时间 n 1500; % 转速 rpm omega 2*pi*n/60; % 角速度 rad/s T_mesh 2*pi/(omega*z1); % 啮合周期 (s) % 时间轴 t linspace(0, 2*T_mesh, 1000); % 双齿叠加第二对齿的刚度曲线是第一对的平移 shift round(eps_a * N); % 平移点数 k_total zeros(size(t)); for j 1:length(t) % 映射到单齿对曲线上的位置 pos mod(t(j)/T_mesh, 1) * N; idx1 max(1, min(N, round(pos)1)); k1 k_single(idx1); % 第二对齿位置偏移一个重合度 pos2 pos eps_a*N; if pos2 N idx2 round(pos2)1; k2 k_single(max(1,min(N,idx2))); k_total(j) 1/(1/k1 1/k2); else k_total(j) k1; end end plot(t*1e3, k_total/1e6, LineWidth, 1.5); xlabel(时间 (ms)); ylabel(啮合刚度 (MN/m)); grid on;逻辑说明单齿对刚度曲线是一条随啮合位置变化的曲线。双齿啮合时第二对齿的啮合位置比第一对超前一个重合度对应的弧长所以把单齿对曲线平移 eps_a*N 个点就是第二对齿的刚度。两对齿并联柔度相加再取倒数。当第二对齿还没进入啮合或已经退出时只有单齿对承载直接取 k1。参数说明shift 是平移点数由重合度决定。eps_a 算出来通常在 1.5 到 1.8 之间对应双齿啮合区占 50% 到 80%。时间轴取两个啮合周期是为了看周期性。刚度单位用 MN/m 是因为钢齿轮单齿对刚度典型值在 1e8 到 1e9 N/m 量级。4. 避坑与排查刚度曲线不对时先查这 5 处4.1 刚度值整体偏大或偏小一个数量级现象跑出来的刚度峰值只有 1e6 N/m 或者高达 1e11 N/m和文献里的 1e8 到 1e9 对不上。原因单位不统一。模数用了毫米但弹性模量用了帕或者齿宽用了毫米但没转成米。势能法里长度单位必须统一到米力单位统一到牛顿。解决在脚本开头把所有几何量统一转成米E 用帕。跑完先看单齿对刚度峰值钢齿轮齿宽 20 mm 时应该在 5e8 N/m 左右。差一个数量级就是单位问题差两个数量级检查积分上下限。4.2 刚度曲线出现周期性毛刺现象曲线整体趋势对但上面叠加了高频锯齿。原因齿廓离散点太少数值积分在齿根附近截面变化剧烈的地方误差大。或者平移叠加时索引取整引入了跳变。解决把 N 从 200 加到 500毛刺通常会消失。如果还有检查平移逻辑里 idx2 的计算用线性插值代替最近点取整。4.3 双齿啮合区刚度反而比单齿区低现象叠加后双齿区的刚度比单齿区还小违反直觉。原因并联公式写反了。两对齿并联是柔度相加再取倒数不是刚度直接相加。如果写成 k1k2双齿区会异常大如果写成 1/(k1k2)双齿区会异常小。解决确认公式是 1/(1/k1 1/k2)。这个式子等价于 k1*k2/(k1k2)永远小于 min(k1,k2)但大于单对齿柔度倒数的一半。4.4 重合度算出来小于 1现象eps_a 算出 0.8 之类导致没有双齿啮合区。原因齿顶圆半径算错或者压力角用错。标准直齿轮重合度一般大于 1.2小于 1 说明齿顶高系数或者压力角输入有误。解决检查 ra r m 是否用了正确的齿顶高系数标准是 1。压力角确认是 20 度不是 25 度。如果确实设计成小重合度那单齿啮合区就是全部叠加逻辑要相应简化。4.5 时间轴和啮合周期对不上现象刚度曲线周期和理论啮合周期不一致或者一个周期里出现了两次波动。原因啮合周期公式用错。啮合周期是转过一个齿距的时间等于 2π/(ω*z)不是 2π/ω。z 用主动轮齿数。解决确认 T_mesh 2pi/(omegaz1)。如果曲线一个周期里有两个波峰说明时间轴取长了或者叠加逻辑里平移量算错导致重复叠加。5. 进阶把刚度曲线喂给动力学模型并验证5.1 从刚度曲线到动力学激励拿到 k_total(t) 只是第一步真正要用它得把它变成动力学方程里的时变系数。常见做法是把刚度曲线做傅里叶展开取前几阶谐波作为激励% 对刚度曲线做 FFT K fft(k_total - mean(k_total)); K_mag abs(K)/length(t)*2; f (0:length(t)-1)/(t(end)-t(1)); % 取前 5 阶谐波 harmonics 1:5; K_harm K_mag(harmonics); fprintf(前5阶刚度谐波幅值 (MN/m):\n); disp(K_harm/1e6);逻辑说明减去均值后做 FFT得到各阶谐波幅值。啮合频率处的幅值就是刚度波动的主要激励。前 5 阶通常够用更高阶幅值衰减很快。这些谐波幅值直接对应动力学方程里刚度矩阵的时变部分。参数说明采样点数要够至少覆盖 10 个啮合周期。如果只取 2 个周期做 FFT频率分辨率太差谐波幅值不准。我一般会跑 20 个周期再截取稳态段做分析。5.2 用能量守恒做一次自检势能法算出来的刚度可以用能量法反推变形量做自检。在啮合点施加一个已知法向力 F用柔度算出总变形 δ F/k然后算变形能 U 0.5Fδ。这个能量应该等于四个势能分量之和。如果对不上说明柔度计算有遗漏。% 自检取齿顶附近一个啮合点 F_test 1000; % N idx_test round(0.8*N); delta_total delta_b(idx_test) delta_s(idx_test) ... delta_a(idx_test) delta_h; k_test 1/delta_total; U_deform 0.5 * F_test^2 * delta_total; % 四个分量能量之和 U_sum 0.5*F_test^2 * (delta_b(idx_test) delta_s(idx_test) ... delta_a(idx_test) delta_h); fprintf(变形能: %.6f J, 分量和: %.6f J\n, U_deform, U_sum);这个自检看起来是恒等式但实际写代码时如果某个柔度分量漏乘了系数或者积分上下限写错这里就会暴露。我习惯每改一次参数就跑一次这个检查比盯着曲线看靠谱。5.3 一个容易被忽略的细节齿根过渡曲线标准渐开线公式只描述工作齿廓齿根到渐开线起点之间有一段过渡曲线。这段曲线对弯曲柔度贡献不小尤其是齿根附近。如果直接用渐开线从齿根圆开始积分会高估齿根刚度。常见做法是用齿条刀具圆角包络出过渡曲线或者用一段圆弧近似。我一般用圆弧近似半径取 0.38m对刚度峰值影响在 5% 以内但曲线形状会更接近实测。跑完这套脚本你会得到一条随啮合位置变化的刚度曲线能直接喂给齿轮动力学模型。参数改起来快试不同齿数、模数、齿宽对刚度波动的影响一个下午能跑几十组。真正花时间的不是写代码是搞清楚每个参数背后的物理意义以及曲线不对时知道去哪找原因。希望帮到你。本文还有配套的精品资源点击获取