简介这份MATLAB脚本面向齿轮动力学建模与振动分析方向的研究生、工程师及科研人员用于求解直齿轮啮合过程中随时间动态变化的啮合刚度。核心算法基于势能法通过建立轮齿变形能与外载荷之间的关系结合齿面接触几何、材料弹性特性及啮合位置变化逐点计算单对齿从啮入到啮出全过程的刚度数值输出时间序列对应的刚度数组可直接用于系统动力学建模、振动响应分析或故障特征提取。资源包共6个文件约42KB以m脚本为主体辅以py脚本、txt说明、png结果图及配置文件结构清晰无需额外工具箱在MATLAB基础平台即可运行。输入参数涵盖模数、齿数、压力角、齿宽、弹性模量、泊松比及转速等基础齿轮参数便于按需调整。目前已有36人学习下载适合希望快速获取时变啮合刚度计算工具、减少重复编程工作的读者参考使用。1. 直齿轮时变啮合刚度到底怎么算一份能跑通的势能法 MATLAB 脚本齿轮箱振动噪声做多了就会发现很多问题的根子不在轴承、不在箱体而在那一对齿面之间随时间跳动的啮合刚度。直齿轮重合度低单齿啮合和双齿啮合交替出现刚度曲线呈明显的方波状跳变这个跳变就是后续动态响应里那些边频带的激励源。想把齿轮动力学仿真做扎实第一步不是去调龙格库塔步长而是先把时变啮合刚度这条曲线算对。这份 MATLAB 脚本用势能法把直齿轮的时变啮合刚度拆成五部分能量分别积分再求和输入只有齿数、模数、压力角、齿宽、弹性模量这几个基本参数输出就是随啮合位置变化的刚度数组。适合正在做齿轮动力学建模、需要自己生成刚度激励的从业者和研究生也适合想搞明白势能法到底在算什么的人。2. 势能法拆解五类能量怎么对应到刚度分量2.1 从单齿悬臂梁到啮合刚度的物理映射势能法的核心思路很朴素把轮齿当成一个变截面的悬臂梁载荷作用在齿面接触点上梁因此发生弯曲、剪切、轴向压缩和基础弹性变形接触点本身还有赫兹接触变形。每一种变形都对应一份弹性势能把五份能量加起来用总能量对法向力求二阶导或者直接用能量与力的关系反推就得到这一对齿在当前位置的啮合刚度。具体拆开看弯曲势能对应齿根弯曲变形剪切势能对应截面剪切变形轴向压缩势能对应齿廓方向的压缩基础弹性势能对应齿根圆角以下那部分基体的变形赫兹接触势能对应两齿面接触区的局部变形。前四项是单齿的柔度贡献第五项是接触对的贡献。单齿柔度是这四项之和一对啮合齿的总柔度是两齿柔度加接触柔度刚度就是柔度的倒数。这里有个容易翻车的地方很多人以为势能法只能算单齿其实它天然适合算一对齿。因为总柔度是可加的两齿各自的弯曲、剪切、轴向、基础柔度加上接触柔度直接得到啮合对的综合柔度。脚本里把主动轮和从动轮分别算柔度再叠加就是这个道理。2.2 脚本的模块划分与数据流拿到脚本先别急着跑花两分钟看清楚它的结构。整个脚本按数据流分成四块参数定义区、单齿柔度计算函数、啮合过程循环、刚度合成与绘图。参数定义区把所有几何和材料参数集中放在开头改参数只动这一块。单齿柔度计算是核心里面按啮合位置逐点积分。啮合过程循环负责判断当前是单齿啮合区还是双齿啮合区决定叠加几个齿对的柔度。最后合成刚度曲线并画图。%% 直齿轮时变啮合刚度计算 - 势能法 % 参数定义区改这里就行 m 3; % 模数 mm z1 20; % 主动轮齿数 z2 40; % 从动轮齿数 alpha 20*pi/180; % 压力角 rad b 20; % 齿宽 mm E 206000; % 弹性模量 MPa nu 0.3; % 泊松比 ha 1; % 齿顶高系数 c 0.25; % 顶隙系数 % 基本几何计算 r1 m*z1/2; r2 m*z2/2; % 分度圆半径 rb1 r1*cos(alpha); rb2 r2*cos(alpha); % 基圆半径 ra1 r1 ha*m; ra2 r2 ha*m; % 齿顶圆半径 rf1 r1 - (hac)*m; rf2 r2 - (hac)*m; % 齿根圆半径这段参数区看着简单但模数和齿数决定了后面所有几何量的量级。模数单位是毫米弹性模量用兆帕最后刚度出来的单位是牛每微米这是齿轮动力学里常用的量纲。如果后面发现刚度值大得离谱或者小得看不见先回来检查这两个单位是不是混了。2.3 单齿柔度的积分实现单齿柔度计算是脚本里最需要看懂的部分。对每一个啮合位置需要沿齿廓从齿根到当前接触点做积分。被积函数里包含截面惯性矩、截面面积、齿廓坐标等随位置变化的量。function [Ub, Us, Ua] single_tooth_energy(y, params) % y: 接触点距齿根的高度 % 返回弯曲、剪切、轴向压缩势能 m params.m; alpha params.alpha; b params.b; E params.E; nu params.nu; G E/(2*(1nu)); % 剪切模量 % 沿齿高离散积分 N 200; h y/N; Ub 0; Us 0; Ua 0; for i 1:N yi i*h; % 截面半角与齿厚 phi alpha - atan(yi/params.r); s 2*params.rb*sin(phi)/cos(alpha); I b*s^3/12; % 截面惯性矩 A b*s; % 截面面积 % 弯曲、剪切、轴向能量密度累加 Ub Ub (yi-yi)^2/(2*E*I)*h; % 简化示意 Us Us 1.2/(2*G*A)*h; Ua Ua 1/(2*E*A)*h; end end上面这段是结构示意实际脚本里被积函数要严格按势能法公式写弯曲项里是力矩的平方除以两倍抗弯刚度剪切项里是剪力平方除以两倍抗剪刚度轴向项类似。参数 N 是积分点数取 200 一般够用齿数多、模数小的时候可以加到 500。G 是剪切模量由弹性模量和泊松比换算。这里最容易出错的是截面半角 phi 的计算它随接触点高度变化算错了整条刚度曲线都会偏。2.4 啮合循环与刚度合成单齿柔度算完接下来要处理啮合过程。直齿轮的重合度一般在 1.2 到 1.8 之间意味着啮合过程中有时是一对齿接触有时是两对齿同时接触。脚本用一个循环遍历主动轮转角每个转角下判断当前有几对齿在啮合把对应的柔度叠加。% 啮合循环 theta linspace(0, 2*pi/z1, 500); % 一个齿距的转角范围 kt zeros(size(theta)); for k 1:length(theta) % 判断当前啮合齿对数 n_pairs judge_contact_pairs(theta(k), params); total_compliance 0; for p 1:n_pairs [Ub1,Us1,Ua1] single_tooth_energy(y1(p), params); [Ub2,Us2,Ua2] single_tooth_energy(y2(p), params); U_total Ub1Us1Ua1Ub2Us2Ua2 U_contact; total_compliance total_compliance 1/(2*U_total); end kt(k) 1/total_compliance; end循环里 judge_contact_pairs 是判断齿对数的关键函数它根据当前转角和重合度决定是一对还是两对。total_compliance 累加的是柔度最后取倒数得到刚度。注意这里每对齿的柔度是两齿柔度加接触柔度接触柔度由赫兹公式给出和法向力有关但势能法里通常取一个等效值。500 个转角采样点对一个齿距来说够密了画出来的刚度曲线足够光滑。3. 参数怎么设模数、齿数、齿宽对刚度曲线的影响3.1 几何参数对平均刚度的影响规律参数不是随便填的每一个都会实打实改变刚度曲线的形状和幅值。模数增大齿厚增加弯曲和剪切柔度下降平均刚度上升。齿数增多在模数不变的情况下分度圆变大单齿更像一根细长梁柔度上升平均刚度下降。齿宽增加直接增加截面面积和惯性矩所有柔度分量都下降刚度近似线性上升。参数变化方向平均刚度变化主要影响分量模数 m增大上升弯曲、剪切齿数 z增大下降弯曲、剪切齿宽 b增大上升全部压力角 alpha增大上升弯曲、轴向弹性模量 E增大上升全部这张表可以当调参的粗略指南。实际做仿真时如果发现刚度幅值比文献里低一个量级先看齿宽是不是填成了毫米但公式里当成了米。如果刚度曲线波动太大检查重合度算得对不对重合度接近整数时单双齿交替会特别剧烈。3.2 重合度与单双齿交替的判断重合度是决定刚度曲线形状的关键。直齿轮重合度等于端面重合度一般在 1.2 到 1.8 之间。重合度小于 1 是不可能的大于 2 意味着始终至少两对齿接触刚度曲线会平缓很多。脚本里判断单双齿交替的逻辑依赖重合度计算。% 重合度计算 eps_alpha (sqrt(ra1^2-rb1^2) sqrt(ra2^2-rb2^2) - (r1r2)*sin(alpha)) / (pi*m*cos(alpha)); fprintf(端面重合度 %.4f\n, eps_alpha); % 单双齿交替判断 if eps_alpha 1 error(重合度小于1齿轮无法连续传动); elseif eps_alpha 2 % 存在单齿啮合区 single_zone_ratio 2 - eps_alpha; else % 始终双齿以上啮合 single_zone_ratio 0; end重合度算出来之后single_zone_ratio 表示一个齿距内单齿啮合所占的比例。这个值越大刚度曲线的方波特征越明显。如果算出来重合度是 1.0 附近单齿区几乎占满整个齿距刚度曲线会呈现很深的谷。这时候要小心实际齿轮加工和安装误差会让这个谷变得不那么尖锐但仿真里如果直接用理想值动态响应会偏大。3.3 材料参数与单位一致性检查材料参数就两个弹性模量和泊松比。钢齿轮弹性模量取 206000 兆帕泊松比 0.3这是默认值。如果换成其他材料比如粉末冶金或者塑料齿轮弹性模量可能低一个量级刚度曲线整体下移。泊松比影响剪切模量对剪切柔度有影响但不如弹性模量敏感。单位一致性是新手最容易翻车的地方。模数用毫米齿宽用毫米弹性模量用兆帕算出来的刚度单位是牛每微米。如果弹性模量误填成帕斯卡刚度会小六个量级。如果齿宽误填成米刚度会大一千倍。脚本开头加一行单位注释每次改参数前看一眼能省很多排查时间。4. 避坑与排查刚度曲线不对时先查这五处4.1 刚度曲线出现负值现象画出来的刚度曲线在某些转角下变成负数。原因柔度叠加时符号搞反了或者接触柔度算成了负值。势能法里所有势能都是正的柔度也是正的刚度取倒数不可能是负的。出现负值一定是中间某一步把能量或者柔度的符号弄错了。解决在循环里加一句判断如果 total_compliance 小于等于零就打印当前转角和相关中间变量。常见的是接触柔度公式里法向力方向取反导致符号错误。把接触柔度单独拿出来验证用一个已知工况手算对比。4.2 刚度幅值比预期低一个量级现象刚度算出来只有几十牛每微米文献里同类齿轮是几百。原因单位不一致。最常见的是弹性模量填成了帕斯卡而不是兆帕或者齿宽填成了米而不是毫米。也有可能是积分点数太少积分没收敛。解决先检查参数区所有量的单位列一个单位对照表。然后把积分点数 N 从 200 加到 1000看刚度值是否收敛。如果加了点数变化不大那就是单位问题。如果变化明显说明原来点数不够。4.3 单双齿交替处曲线跳变过于剧烈现象刚度曲线在单双齿切换点出现垂直跳变像台阶一样。原因理想模型下确实存在跳变但实际齿轮由于齿廓修形和安装误差跳变会被平滑掉。如果仿真里直接用理想跳变后续动力学响应会偏大。解决可以在切换点附近加一个过渡区用平滑函数把两段刚度曲线过渡一下。常见做法是用一个宽度为几个采样点的余弦过渡。这不是势能法本身的要求是为了让仿真更接近实际。4.4 重合度算出来大于 2 或小于 1现象脚本报错或者刚度曲线形状明显不对。原因齿顶圆半径或者基圆半径算错了。齿顶圆半径等于分度圆半径加齿顶高齿顶高等于齿顶高系数乘模数。基圆半径等于分度圆半径乘压力角余弦。这两个算错重合度就错了。解决把齿顶圆半径、基圆半径、分度圆半径打印出来手算核对。标准直齿轮齿顶高系数取 1顶隙系数取 0.25。如果用了变位齿轮齿顶圆半径公式要改脚本里默认是标准齿轮。4.5 积分不收敛或结果随积分点数剧烈变化现象把积分点数从 200 改到 500刚度值变化超过百分之十。原因被积函数在齿根附近变化剧烈固定步长积分不够。齿根处截面惯性矩很小被积函数有尖峰。解决把等步长积分改成自适应积分或者在齿根附近加密采样。简单做法是把积分点数直接加到 1000 以上代价是计算时间增加。如果只是做几次仿真多花几秒无所谓。如果要嵌入优化循环就得考虑用解析积分或者查表法加速。5. 从刚度曲线到动力学激励一个验证技巧刚度曲线算出来只是第一步真正要用它得把它变成动力学模型里的激励。这里分享一个验证技巧把算出来的时变啮合刚度做傅里叶变换看它的谐波成分。直齿轮的时变啮合刚度主要包含啮合频率及其倍频如果傅里叶谱里出现了非啮合频率的杂波说明刚度曲线在单双齿切换处有数值毛刺需要平滑。% 对刚度曲线做频谱分析验证 kt_mean mean(kt); kt_ac kt - kt_mean; % 去掉直流分量 Fs length(theta); % 采样点数 Y fft(kt_ac); P2 abs(Y/Fs); P1 P2(1:Fs/21); P1(2:end-1) 2*P1(2:end-1); f (0:(Fs/2))/Fs; plot(f, P1); xlabel(归一化频率); ylabel(幅值); title(时变啮合刚度频谱);这段代码把刚度曲线去掉均值后做快速傅里叶变换看幅值谱。正常情况下谱线应该集中在啮合频率的整数倍上。如果看到很多分散的谱线说明刚度曲线有数值噪声。这时候回去检查积分点数和单双齿切换处的处理。我自己的习惯是每次改完参数先看刚度曲线形状再看频谱两个都对了才拿去做动力学仿真。有一次偷懒没看频谱直接跑动力学结果响应里多了一堆莫名其妙的频率成分查了半天才发现是刚度曲线在切换点有毛刺。从那以后我每次生成刚度曲线都强制走一遍频谱检查花不了一分钟能省掉后面几个小时的排查。希望帮到你。本文还有配套的精品资源点击获取