在齿轮箱的振动噪声仿真里时变啮合刚度TVMS绕不开。直齿轮的刚度计算有很多现成文献但换成斜齿轮事情就复杂了接触线是斜着扫过整个齿宽的单对齿和多对齿交替承载的规律跟直齿完全不是一个逻辑。这时候要是不管三七二十一拿直齿轮的算法硬套斜齿轮算出来的刚度曲线不光幅值不对趋势都是错的。所以我这些年做齿轮传动仿真一直是用切片法把所有斜齿轮问题降维成直齿轮问题再配合势能法挨个算每一片薄齿的啮合刚度最后沿接触线累加出整条时变刚度曲线。这篇就把这套完整的Matlab实现思路、关键公式和调试经验一次说透。这篇内容比较适合正在做齿轮动力学、箱体振动噪声、传动系统NVH的同行也适合研究生刚接触齿轮仿真、被各种啮合刚度文献绕晕的新手。只要你装了MatlabR2019b以后的版本都够用跟着下面的几何推导和代码结构走是能直接复现出一条靠谱的斜齿轮时变啮合刚度曲线的。1. 先搞明白为什么斜齿轮必须用切片法来算刚度1.1 直齿和斜齿的核心差别直齿轮啮合时一对齿的接触线是一条与轴线平行的直线整个齿宽在同一个瞬间同时进入啮合。单齿啮合区、双齿啮合区划分得非常清楚刚度跳变也很明显。斜齿轮不一样。螺旋角让齿面错开了接触线从齿顶的一个角开始斜着往齿根方向扫过去接触线长度是先增后减的而且每一瞬间参与啮合的轮齿对数也在连续变化。这导致斜齿轮的时变啮合刚度波动明显比直齿轮小但绝不是没有波动只是波动形态更复杂跟螺旋角、齿宽、重合度都强相关。这里有个关键概念斜齿轮的总重合度是端面重合度加轴向重合度也就是 εγ εα εβ。其中轴向重合度 εβ b·sinβ/(π·m_n)b 是齿宽β 是螺旋角m_n 是法面模数。轴向重合度小于1的时候接触线是断续的会有啮合空档大于1的时候接触线连续覆盖刚度曲线更平滑。做切片法之前先把这个值算清楚否则后面相位对不齐结果就是一团糟。1.2 切片法的物理图像把斜齿轮切成一堆直齿轮切片法说白了就是微积分思想的工程化应用。把斜齿轮沿轴线切成 N 个等厚度的薄圆盘每个圆盘可以近似看成一个极薄的直齿轮。由于螺旋角存在每一个薄片在圆周方向相对于相邻薄片都错开了一个微小的角度也就是相位差。切片越薄这个直齿轮近似就越准。于是斜齿轮的啮合刚度问题就转化成一组“错开相位、厚度相同”的直齿轮薄片的刚度并联问题。在某一时刻 t每个薄片对应不同的啮合位置对接触线都有局部刚度贡献把所有同时参与啮合的薄片刚度累加就能得到整个斜齿轮副在该时刻的总啮合刚度。这个累加过程本质就是在沿接触线积分所以切片法天然适合斜齿轮这种接触线倾斜的几何。你可能会问为什么不用三维有限元直接算能算但是慢。齿轮啮合刚度分析要扫一个完整啮合周期通常分 30 到 60 步每步都是一次非线性接触分析计算成本非常可观。切片法的每个片是纯解析公式Matlab跑几百个切片也是秒级完成做参数扫描优化特别合适。1.3 势能法为什么和切片法是天生一对势能法也叫能量法把轮齿看成悬臂梁认为啮合力做的功转化为五部分变形能弯曲势能、剪切势能、轴向压缩势能、赫兹接触势能、齿基体弹性变形势能。每一部分都可以用一个积分公式表示最后把五个刚度分量并联得到单齿啮合刚度。它和切片法配合起来非常自然。每个薄片都是一个直齿轮模型直齿轮的势能法公式在二维截面里定义得好好的切片厚度方向的泊松约束近似不影响公式主项。你只需要把每个切片当独立的悬臂梁算一遍再把结果按接触顺序叠加。而且势能法对渐开线齿廓和齿根过渡曲线的建模误差在切片厚度足够薄、齿廓离散点足够密时可以压到很小经验上跟精细有限元解误差大概率能控制在5%以内。2. 五种势能分量的公式与几何量准备2.1 从法面参数到端面参数先解决单位制问题斜齿轮设计图纸给的通常是法面模数 m_n、法面压力角 α_n、螺旋角 β以及齿数 z。但在势能法公式里所有几何量都必须在端面坐标系里算不然齿廓坐标和宽度方向就串了。端面模数的换算关系很简单m_t m_n / cosβ端面压力角 tan α_t tan α_n / cosβ分度圆直径 d m_t · z基圆直径 d_b d · cosα_t端面齿顶高系数 h_a* 和顶隙系数 c* 也要按端面方向投影修正但多数文献里仍用法面值近似我踩过最狠的一个坑就是角度量纲。Matlab里 sin、cos 默认输入是弧度你从Excel抄个压力角20直接丢进公式结果荒谬到天际。高端齿轮设计软件里压力角单位设置不一样自己写程序时最好统一规定输入参数用度进入计算函数第一行就全部转弧度内层只用弧度运算输出角度再转回度。2.2 渐开线齿廓与齿根几何的建模势能法积分的核心是把齿廓离散成很多个点每个点对应一个截面位置然后逐点计算截面惯性矩和截面积。所以先得生成渐开线齿廓。渐开线的极坐标参数方程如下基圆半径 r_b压力角变量 θ其实就是展开角极径 r r_b / cosθ极角 φ tanθ - θ这是渐开线函数 invθ在Matlab里离散时可以这样写从齿根圆半径 r_f 到齿顶圆半径 r_a 等间距取 θ然后转成笛卡尔坐标。齿根过渡曲线那一段严格说是圆弧加直线但在切片法里过渡曲线对弯曲变形的影响远小于渐开线段很多工程文章直接把它简化成一段圆弧圆弧半径取 0.38·m_n 左右影响可接受。我建议齿廓采样点不少于 200 个。太少的时候靠近齿根附近的截面惯性矩会出现锯齿状跳变算出来的弯曲刚度跟着抖太多则慢但也没必要超过 1000毕竟后面还有切片维度总的计算量是两者相乘。2.3 弯曲、剪切、压缩、赫兹、基体五部分刚度公式现在进入正题把五种势能的刚度公式全部列出来。假设啮合力 F 作用在啮合点上啮合点距离齿根截面危险截面的水平距离是 x整个齿根悬臂长度是 d。在齿廓离散位置 x_i 处截面厚度 h_x 可以从几何关系直接求出截面积 A_x 和惯性矩 I_x 就很直观了A_x 2·h_x·LL是切片厚度对于斜齿轮就是 dzI_x L·(2·h_x)^3 / 12弯曲势能刚度k_b 1 / ∫₀ᵈ [x²/(E·I_x)] dx剪切势能刚度k_s 1 / ∫₀ᵈ [1.2/(G·A_x)] dx其中 1.2 是矩形截面的形状系数G 是剪切模量G E/[2(1ν)]。轴向压缩势能刚度k_a 1 / ∫₀ᵈ [sin²α/(E·A_x)] dx这个 α 是啮合力与齿宽端面的夹角通常取端面压力角近似。赫兹接触刚度k_h π·E·L / [4(1-ν²)]齿基体弹性变形刚度用 Sainsot 公式k_f 1 / { L·cos²α / E · [ (L*)²·(u_f/S_f) M*·(u_f/S_f) P*·(1Q*·tan²α) ] }这里的 u_f 是啮合点到齿根圆角中点的距离S_f 是齿根圆弧长度L*、M*、P*、Q* 是一组查表系数它们本身跟齿根圆弧与齿厚比值相关。我常用的系数组如下这组数据来自工业标准文献可以直接查表用系数数值第一组适用常见标准齿L*0.5872M*0.4254P*1.0793Q*0.8330切片厚度 L 在每片里其实就是 dz按我 3.2 节的方式传入不要搞混。五个刚度分量是串联关系所以单齿对的瞬时啮合刚度1/k_total 1/k_h 1/k_b 1/k_s 1/k_a 1/k_f这里的几何量必须对应同一个切片、同一个时刻。啮合点在齿面上移动x、d、u_f 全在变所以单齿刚度是关于啮合位置的函数。这也是“时变”的根本来源。3. Matlab程序架构与关键代码3.1 程序流程从参数定义到刚度曲线整个程序我习惯分成四层。第一层是输入参数定义包括材料参数、齿形几何参数、螺旋角和切片数。第二层是几何预处理计算端面参数、生成渐开线齿廓点、离散断面。第三层是核心仿真循环按时间步进每个时刻扫描所有参与啮合的切片每个切片单独算单齿耦合刚度再累加。第四层是后处理画刚度曲线、做傅里叶谐波分析、跟有限元结果对比。主程序骨架如下注释写清楚每一步的意图%% 参数定义 E 2.06e5; % 弹性模量 MPa nu 0.3; % 泊松比 mn 2.5; % 法面模数 mm zn 23; % 主动轮齿数 zg 57; % 从动轮齿数 beta 24; % 螺旋角 deg alpha_n 20; % 法面压力角 deg b 40; % 齿宽 mm slice_n round(b / 0.5); % 切片数每片0.5mm rot_deg_step 0.5; % 转角步长 deg %% 端面几何换算 mt mn / cosd(beta); alpha_t atand(tand(alpha_n) / cosd(beta)); d1 mt * zn; d2 mt * zg; db1 d1 * cosd(alpha_t); db2 d2 * cosd(alpha_t); ra1 d1/2 mn; % 齿顶圆半径简单近似 rf1 d1/2 - 1.25*mn; % 齿根圆半径近似 %% 几何预处理 [profile_x, profile_h] tooth_profile(mt, alpha_t, z, ...); % 得到每个截面位置和半齿厚 %% 时变刚度主循环 mesh_angle 0 : rot_deg_step : 360/zg; % 一个啮合周期 for k 1 : length(mesh_angle) k_total(k) assemble_slices(mesh_angle(k), ...); end %% 后处理 plot(mesh_angle, k_total);注意一个啮合周期的角度范围严格说是从 360/z 除以总重合度也就是主动轮转过一个基节对应的时间。上面骨架里写 360/zg 是给从动轮的简化写法实际程序中要按啮合线位置计算更精确的起止点后面 3.3 节会说明。3.2 单齿切片刚度函数核心中的核心每个切片的刚度计算我封装成一个独立的函数。它的输入是切片当前对应的啮合位置输出是一个数值该切片在这一瞬间的有效啮合刚度。function k_slice slice_dynamic_stiffness(geo, mesh_x) % geo: 结构体, 包含齿廓离散坐标、材料参数 % mesh_x: 当前啮合点沿啮合线方向的相对位置, mm % 1. 根据mesh_x确定啮合点在齿面上的位置, 更新x,d,u_f % 2. 沿着齿廓积分算五个刚度分量 x geo.x_nodes; % 截面距啮合点的水平距离 hx geo.h_nodes/2; % 半齿厚, mm Ls geo.slice_width; % 切片厚度 mm E geo.E; G geo.G; nu geo.nu; Ax 2 * hx * Ls; % 截面积 Ix Ls .* (2*hx).^3 / 12; % 惯性矩 inv_kb trapz(x, x.^2 ./ (E*Ix)); % 弯曲柔度 inv_ks trapz(x, 1.2 ./ (G*Ax)); % 剪切柔度 inv_ka trapz(x, sin(geo.alpha_t)^2 ./ (E*Ax)); % 轴向压缩 kb 1/inv_kb; ks 1/inv_ks; ka 1/inv_ka; kh pi*E*Ls / (4*(1-nu^2)); % 赫兹刚度 % 基体刚度用Sainsot公式, 这里简化计算u_f和S_f u_f geo.uf_table(mesh_x); S_f geo.Sf; Lstar 0.5872; Mstar 0.4254; Pstar 1.0793; Qstar 0.8330; coef Lstar*(u_f/S_f)^2 Mstar*(u_f/S_f) Pstar*(1Qstar*tan(geo.alpha_t)^2); kf E*Ls / (coef * cos(geo.alpha_t)^2); % 五个分量串联 k_slice 1/(1/kh 1/kb 1/ks 1/ka 1/kf); end这个函数里最需要注意的是积分变量 x 的物理含义。很多文献里弯曲柔度积分是从危险截面到啮合点也就是齿根到受力点方向与齿高方向一致。我在代码里把 x 数组定义成从齿根零点到啮合点的离散距离这样 trapz 直接用不容易出错。人生的教训是不要为了把公式抄齐BSI文献里的符号反而把自己的索引搞乱。3.3 啮合时序与刚度叠加相位错开是灵魂切片法最关键的步骤就是给每一个切片刻画相位差。斜齿轮总有螺旋角 β所以轴向坐标 z 处的切片相对齿轮端面起始位置要额外转过一个角度Δφ(z) (z·tanβ) / (d/2) · 180/π这里 z 是该切片中心相对齿宽一端的轴向距离。也就是说最左端切片的相位差是0最右端切片的相位差最大。要正确用时变刚度必须把这个相位差叠加到该切片自己的“啮合进程角”上。这样一来主循环里的时间步就会对应一幅“切片群像”有的切片刚进入啮合刚度从小变大有的切片在单齿区刚度稳定在最大值还有的切片快要退出刚度衰减。把所有切片同一时刻的 k_slice 加起来就得到整副齿轮副在该转角位置的总啮合刚度。下面是装配函数的核心逻辑function k_total assemble_slices(mesh_rot, geo, phase_offsets) k_total 0; for i 1 : geo.slice_n local_rot mesh_rot phase_offsets(i); % 相位补偿 if local_rot geo.entry_angle local_rot geo.exit_angle mesh_x geo.mesh_x_map(local_rot); % 转角 - 啮合线位置 k_total k_total slice_dynamic_stiffness(geo, mesh_x); else % 该切片当前不参与啮合 continue; end end endentry_angle 和 exit_angle 是单对轮齿在啮合线上的进入和退出角度可以由重合度推算。如果你把 phase_offsets 均匀从 0 分布到总相位差再把每个时刻所有切片遍历一遍就相当于数值实现了沿接触线的积分。这里有个小技巧不要在每个时刻从零开始重新计算所有切片的 stiffness性能建议是先算好一张“切片位置-单齿刚度”的插值表之后按转角查表累加。对参数扫描来说这个优化能省好几倍时间数据量大时非常明显。3.4 结果可视化与傅里叶分析主循环跑完得到一条完整的刚度-转角曲线。接下来要干的活有两个一是画图检查曲线形态二是做FFT分析谐波成分因为齿轮动力学里啮合刚度的谐波直接决定振动激励强度。画图简单plot(转角,刚度) 就行注意横坐标统一成齿轮转角或者啮合周期百分比。我更推荐用“啮合周期比”作为横轴0到1代表完整啮合周期这样方便和文献对比。FFT分析点要注意采样步长。如果你的转角步长是 0.2°齿轮一周有1800个采样点FFT的前几阶谐波都能量得出来。示例代码如下fs 360 / mesh_step; % 每转采样频率(采样点数/转) K_fft fft(k_total); f_axis (0:length(K_fft)-1) * fs / length(K_fft); stem(f_axis(1:20), abs(K_fft(1:20))); xlabel(谐波阶次); ylabel(频谱幅值); title(啮合刚度曲线谐波分析);我一般会把一阶谐波幅值和平均值标出来这两个值在后续齿轮啮合噪声预测里就是要用的激励源参数。如果算出来的刚度曲线平均数值差得太离谱多半是几何层出了问题先别急着怀疑程序回头查齿廓离散和单位。4. 工程上的坑与调试经验4.1 常见的六类错误快查表下面这六类错误我基本都踩过随便列出来你看一眼可能就帮你省一个通宵现象原因排查方法刚度曲线幅值整体偏小几个量级单位不统一混用了米和毫米检查E是否为MPa量级长度是否全为mm曲线出现周期性锯齿且与切片数相同切片数太少阶梯效应增加切片数到每毫米至少2片曲线严重不对称左右端波形不同相位偏移公式中用了分度圆直径而非基圆直径检查Δφ公式里的半径用错了没单齿啮合区消失曲线始终在叠加重合度计算错误边界搞错先用公式算εα和εβ核对刚度突变位置与理论啮合周期不符进入/退出角度设定错误打印调试变量检查entry/exit angle结果和文献对不上但趋势对齿根过渡曲线简化过于粗糙把简化改成标准圆弧加直线还有保留上面表格里每一项我用图文标注过不止一次所以奉劝大家出了问题先跑这套排查清单不要在整代码层面反复怀疑人生。4.2 切片数怎么选收敛性和锯齿效应切片数是个必须调的超参数少了激振锯齿多多了算得慢没有绝对标准。我的经验法则是切片厚度不超过 0.5mm。比如齿宽 40mm就取80片上下如果齿宽很小比如 10mm也可以最少20片。我自己做过一次收敛性试验固定斜齿轮参数不变把切片厚度分别取 2mm、1mm、0.5mm、0.25mm对比刚度平均值和峰峰值切片厚度 mm平均刚度 N/μm峰峰值波动 %计算耗时2.018.36.20.2s1.018.74.50.5s0.518.93.81.1s0.2518.93.72.2s从这个结果能看出0.5mm之后平均刚度基本收敛峰峰值变化也不大了。所以日常计算我固定用 0.5mm 切片既不会太慢也足够描述斜齿轮的人字齿效应。如果你要跟FEM对标再加一个 0.25mm 工况验证收敛即可。4.3 推荐的验证与对标手段别急着把切片法结果当成真值。我通常会做三层次验证第一是解析公式粗校验把斜齿轮螺旋角清零退化回直齿轮再跟 Ma 和 Li 等人发表的直齿解析解对比如果对不上说明基线程序有虫子。第二是有限元静态对标随便用 ANSYS Workbench 或者 ABAQUS 做一个单齿对接触模型施加单位法向力测综合变形跟程序比刚度。第三是利用重合度做动态校验仿真出来的刚度曲线峰值谷值出现的位置必须和理论重合度推算的啮合区间一致。有限元对标的网格密度也要注意齿轮接触区网格尺寸至少要到 0.1mm 以下否则赫兹接触刚度会被网格过大严重低估。我没有说一定要用非线性接触算法线性静力分析配合赫兹修正往往也是一种足够好的对标方式。做这层验证还有一个额外好处把切片法程序和有限元结果一起打包存档写论文或者做报告时直接引用对标图比空口说“程序已验证”有说服力得多。说点个人体会。第一次写这个程序时我以为最难的是势能法公式推导结果真正难的是把几何离散和切片相位理清楚。Matlab天然适合这类矩阵化计算但如果你一开始就用向量化思维把每一片每一时刻的刚度直接组织成二维矩阵来算代码会清爽很多。后来我还把这个程序包成了参数扫描版本输入一系列螺旋角和齿宽一键输出刚度曲线族和对应谐波能量做齿轮宏参优化非常高效。你照着这个结构写后面扩展成双齿啮合顺序分析、考虑齿廓修形都会省事很多。