企业官网建设流程全解析

简介：这份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 - (ha+c)*m; rf2 = r2 - (ha+c)*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*(1+nu)); % 剪切模量 % 沿齿高离散积分 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 = Ub1+Us1+Ua1+Ub2+Us2+Ua2 + 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 几何参数对平均刚度的影响规律

参数不是随便填的，每一个都会实打实改变刚度曲线的形状和幅值。模数增大，齿厚增加，弯曲和剪切柔度下降，平均刚度上升。齿数增多，在模数不变的情况下分度圆变大，单齿更像一根细长梁，柔度上升，平均刚度下降。齿宽增加直接增加截面面积和惯性矩，所有柔度分量都下降，刚度近似线性上升。

这张表可以当调参的粗略指南。实际做仿真时，如果发现刚度幅值比文献里低一个量级，先看齿宽是不是填成了毫米但公式里当成了米。如果刚度曲线波动太大，检查重合度算得对不对，重合度接近整数时单双齿交替会特别剧烈。

3.2 重合度与单双齿交替的判断

重合度是决定刚度曲线形状的关键。直齿轮重合度等于端面重合度，一般在 1.2 到 1.8 之间。重合度小于 1 是不可能的，大于 2 意味着始终至少两对齿接触，刚度曲线会平缓很多。脚本里判断单双齿交替的逻辑依赖重合度计算。

% 重合度计算 eps_alpha = (sqrt(ra1^2-rb1^2) + sqrt(ra2^2-rb2^2) - (r1+r2)*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/2+1); P1(2:end-1) = 2*P1(2:end-1); f = (0:(Fs/2))/Fs; plot(f, P1); xlabel('归一化频率'); ylabel('幅值'); title('时变啮合刚度频谱');

这段代码把刚度曲线去掉均值后做快速傅里叶变换，看幅值谱。正常情况下，谱线应该集中在啮合频率的整数倍上。如果看到很多分散的谱线，说明刚度曲线有数值噪声。这时候回去检查积分点数和单双齿切换处的处理。

我自己的习惯是每次改完参数，先看刚度曲线形状，再看频谱，两个都对了才拿去做动力学仿真。有一次偷懒没看频谱，直接跑动力学，结果响应里多了一堆莫名其妙的频率成分，查了半天才发现是刚度曲线在切换点有毛刺。从那以后我每次生成刚度曲线都强制走一遍频谱检查，花不了一分钟，能省掉后面几个小时的排查。希望帮到你。

本文还有配套的精品资源，点击获取