企业官网建设流程全解析

简介：面向机械工程与齿轮动力学方向的 MATLAB 源码包，围绕单自由度直齿轮副非线性动力学问题，提供从动力学方程建立、数值求解到结果可视化的完整代码实现。资源共 1 个 m 文件，压缩包大小仅 1KB，虽为精简脚本，但涵盖相图绘制与傅里叶变换分析等关键环节，适合学习齿轮系统振动特性、开展课程设计或科研预研的读者参考。目前已有 678 人学习/下载。通过阅读与运行这套代码，可掌握单自由度齿轮系统建模思路，理解啮合冲击、弹性变形等因素对动态响应的影响，并借助 MATLAB 实现相空间轨迹与频域特征提取，为后续更复杂的多自由度齿轮动力学研究打下基础。

1. 为什么拆这个单自由度齿轮动力学MATLAB程序？

拿到AnalysisofNolinearDynamicsinaSpurGearPairSystem.m这个文件，我第一反应以为又是常规的齿轮箱振动分析。真正跑起来之后才发现，它把直齿轮副的啮合过程压缩成一个单自由度模型，只保留啮合线方向的相对位移，然后用相图和傅里叶变换两条路径去观察系统状态。这个思路在工程上很有用：建模过程不堆自由度，物理意义直观；求解速度快，能反复改参数看趋势；非线性行为如齿侧间隙引起的冲击和脱啮，也能在这个低维模型里体现出来。无论你是在做齿轮传动设计、状态监测还是故障诊断，这套从方程到频域分析的链路都值得拆开读一遍。它不追求和有限元结果逐点一致，但能快速回答“参数改变后系统会稳定还是分岔”这类关键问题。

2. 单自由度齿轮模型的理论基础：从牛顿第二定律到啮合方程

2.1 为什么齿轮副能简化成单自由度

直齿轮副啮合时，主要振动集中在啮合线方向。如果把两个齿轮和支承轴等效为集中质量，啮合区域等效为刚度和阻尼随时间变化的弹簧，那么系统的运动就可以用沿啮合线方向的相对位移一个变量来描述。这就是单自由度动力学方程（SDOF）的由来。

这种简化不是为了偷懒。在实际齿轮箱中，多自由度系统里各模态相互耦合，参数标定困难，很多时候连模态振型都得不到一致结论，而单自由度模型能抓住最关键的啮合刚度激励。比如重合度在1~2之间时，啮合齿对数在1和2之间变化，啮合刚度周期性波动，这个激励源决定了振动的基本频率成分。先把这个频率成分分析清楚，再考虑轴向扭转等因素，是工程上成熟的递进思路。

2.2 齿轮动力学方程与关键参数

基于牛顿第二定律，直齿轮副啮合线方向的动力学方程可以写成：

m \ddot{x} + c \dot{x} + k(t) f(x) = F_m + F_h cos(ω t)

其中，x 是啮合线方向的相对位移，m 是等效质量，c 是啮合阻尼，k(t) 是时变啮合刚度（简化为平均刚度加一次谐波），f(x) 表示齿侧间隙引起的恢复力非线性，F_m 是平均载荷，F_h cos(ωt) 是啮合冲击激励。ω 即啮合频率，等于齿轮转频乘以齿数。

这里把方程参数整理成一张表，方便写代码时对照：

注意，f(x) 是典型的非光滑函数：

f(x) = x - b, 当 x > b；0, 当 |x| ≤ b；x + b, 当 x < -b。

这是单自由度模型中非线性最主要的来源。当振动位移小于侧隙时，齿轮处于脱啮状态，此时两个齿面没有接触，恢复力为零，系统的刚度瞬时下降，产生冲击。这种非线性在总装传动中会引出一系列周期倍化和混沌现象。

2.3 非线性从哪里来：时变刚度与冲击的耦合

齿轮动力学中的非线性不只是齿侧间隙。时变啮合刚度本身是周期性时变参数，当振动幅度足够大，齿面脱离接触时，间隙函数开始起作用，两个非线性机制叠加在一起，系统就可能在某一转速区间出现跳跃、次谐波共振甚至混沌。

在MATLAB求解时，最方便的是把这种非光滑函数写成条件判断。每一个时间步都重新判断x的位置，从而决定恢复力的表达式。由于刚度突变，数值积分需要使用对刚性不敏感的求解器，常见做法是先用ode45试算，如果在间隙边界附近出现收敛慢或振荡，就换用ode15s。源码包中的脚本主要用ode45，因为齿轮系统的质量、阻尼、刚度参数通常在非刚性范围内。

这一章我们建立的理论模型是整个分析的基石。下一章直接看代码，看方程如何变成可执行的MATLAB脚本。

3. MATLAB求解：把方程写进脚本，跑出位移和速度

3.1 二阶方程降阶为一阶ODE方程组

MATLAB的ode45只能处理一阶常微分方程组，所以先把二阶方程改写成状态空间形式。令 z1 = x，z2 = dx/dt，则得到：

dz1/dt = z2
dz2/dt = (F_m + F_h cos(ωt) - c z2 - k(t) f(z1)) / m

这个形式直接对应函数文件的输入输出。先看一下齿轮微分方程函数怎么写：

function dz = gearODE(t, z, param) % 单自由度直齿轮动力学方程 % z(1) = x（啮合线方向位移） % z(2) = dx/dt（速度） % param 为结构体，包含全部系统参数 m = param.m; % 等效质量 kg c = param.c; % 啮合阻尼 N.s/m km = param.km; % 平均啮合刚度 N/m ka = param.ka; % 刚度波动幅值 N/m w = param.w; % 啮合频率 rad/s b = param.b; % 齿侧间隙半宽 m Fm = param.Fm; % 平均载荷 N Fh = param.Fh; % 冲击激励幅值 N % 时变啮合刚度：平均刚度 + 一次谐波 kt = km + ka * cos(w * t); % 间隙非线性函数 if z(1) > b gap = z(1) - b; elseif z(1) < -b gap = z(1) + b; else gap = 0; end % 状态方程 dz = zeros(2,1); dz(1) = z(2); dz(2) = (Fm + Fh * cos(w * t) - c * z(2) - kt * gap) / m; end

注意这段代码中，kt在每一步重新计算，因为它是时间t的函数。gap根据当前位移判断是否处于啮合状态。当位移落在间隙区内时，恢复力设为零，这就实现了脱啮过程的模拟。如果直接把gap写成一个符号分段函数再调用subs，数值积分每一步都会做符号计算，速度慢一个数量级，我一般在正式脚本里直接用if-else。

3.2 主脚本中如何调用ode45

定义好方程函数后，主脚本里先给出系统参数。以一组典型直齿轮参数为例：

% 系统参数 param.m = 2.0; % 等效质量 kg param.c = 80; % 阻尼系数 N.s/m param.km = 2e6; % 平均啮合刚度 N/m param.ka = 0.6e6; % 刚度波动幅值 N/m param.w = 1885; % 啮合频率 rad/s（转频300Hz × 齿数40 / 2π? 这里直接给角频率） param.b = 40e-6; % 齿侧间隙半宽 m param.Fm = 2500; % 平均载荷 N param.Fh = 500; % 动态载荷幅值 N % 初始条件：位移稍微偏离平衡，速度初值为0 z0 = [1e-5; 0]; % 时间范围：让系统经过瞬态进入稳态 tspan = [0 0.5]; % 求解 [t, z] = ode45(@(t, z) gearODE(t, z, param), tspan, z0);

代码里把啮合频率直接定义成了角频率param.w，方便和方程中的cos(w*t)对应。如果你习惯用转速和齿数计算，可以从w = 2*pi* n_rpm/60 * z得到，这里为突出方程本身而省略换算过程。

时间范围[0 0.5]表示积分0.5秒。假设啮合频率300Hz，0.5秒内有150个啮合周期，足够让瞬态衰减并观察稳态轨迹。如果只关心周期响应，可适当压缩时间；如果做分岔分析，则需要跳过瞬态，只取后半段数据。

ode45默认自动调整步长，返回的t和z是不等间隔的数据点。这里有个常见坑：后面对信号做FFT时需要等间隔采样，所以应从t中重新构建采样频率，或者用interp1重采样。我一般先直接看t(2)-t(1)，确认平均步长是否满足采样要求。

3.3 画相图：观察周期解和混沌

相图把速度作为纵轴、位移作为横轴，把系统的状态轨迹画在同一张图上。它对周期运动和非周期运动的分辨力很强：周期1运动对应一条闭合曲线，周期2是双环闭合，混沌则是填充一定区域的不规则轨迹。

% 画相图 figure('Color', 'w'); plot(z(:,1)*1e6, z(:,2), 'b.-', 'MarkerSize', 2); xlabel('位移 x (μm)'); ylabel('速度 dx/dt (m/s)'); title('单自由度齿轮系统相图'); grid on;

这里把位移从米换算成微米显示，更符合工程读数习惯。坐标轴不要自动缩放过头，否则周期运动的轨迹会贴成一条粗线。建议先把前20%的瞬态数据去掉，只画稳态段：

steady_start = round(0.3 * length(t)); plot(z(steady_start:end,1)*1e6, z(steady_start:end,2), 'b.');

通过试算不同的阻尼和间隙值，你会看到相图从一条细椭圆逐渐变成带毛刺的环，最后变成一大团点云，这就是系统从周期到混沌的演化过程。很多论文里常用的“单自由度系统分岔图”，本质就是把这个过程按某个参数连续扫描再叠加。

4. 傅里叶变换：从时域信号提取啮合频率特征

4.1 为什么用FFT而不是直接看时域

齿轮振动信号里包含周期性的啮合冲击，时域波形只能看出大致的振幅变化，但分不清具体的频率成分。比如一个齿面发生局部剥落，时域上只在对应转角处出现一个尖峰，频域上则表现为啮合频率两侧出现以转频为间隔的边带。为了提取这些特征，必须把信号变换到频域。MATLAB里最常用的就是fft函数。

如果手里已有从实验测得的振动数据，常见做法是先导入CSV文件：

data = readmatrix('vibration_data.csv'); % 第一列时间，第二列加速度 t_data = data(:,1); acc = data(:,2);

然后计算采样频率fs = 1 / (t_data(2) - t_data(1))。但要注意，实际采集数据的采样频率可能不是恒定的，最好先用diff(t_data)看时间间隔有无抖动，有抖动就先重采样。

4.2 用fft分析仿真得到的位移信号

我们直接对上一章求解出的位移信号做FFT。由于ode45返回的是非等间隔时间轴，需要先重采样。常见做法是用t构造新的等间隔时间向量：

fs = 1 / (t(2) - t(1)); % 近似采样频率，ode45步长不均匀 t_uniform = linspace(0, t(end), length(t)); x_uniform = interp1(t, z(:,1), t_uniform, 'linear'); N = length(x_uniform); Y = fft(x_uniform); P2 = abs(Y / N); P1 = P2(1:floor(N/2)+1); P1(2:end-1) = 2 * P1(2:end-1); f = fs * (0:floor(N/2)) / N; figure('Color','w'); plot(f, P1); xlabel('频率 (Hz)'); ylabel('幅值'); xlim([0 1500]); % 只显示0~1500Hz grid on;

代码解释：interp1把非等间隔信号映射到等间隔时间轴上，这一步不能省，否则FFT的频谱会出现虚假频率。P2 = abs(Y/N)做幅值归一化，因为MATLAB的FFT结果数值与信号长度成正比。P1(2:end-1) = 2*P1(2:end-1)是单边谱修正，把负频率的能量折回正频。最后频点序列f从0到奈奎斯特频率，即采样频率的一半。

在实际操作中，如果原始采样频率很高，建议先对信号做带通滤波再FFT，否则低频分量会压过感兴趣的啮合频率。MATLAB里可以用bandpass(x, [500 2000], fs)，简单直接。

4.3 边带分析：齿轮局部故障的指纹

识别齿轮状态有个关键技巧：观察啮合频率附近有没有边带。啮合频率 f_z = z × f_r，其中 z 为齿数，f_r 为转频。边带出现在 f_z ± n·f_r 处，n一般取1到3，这些边带是由齿面缺陷引起的幅值调制或频率调制造成的。

通过对比仿真信号的FFT结果，你会发现当齿侧间隙从20μm增大到80μm时，频谱中f_z附近出现大量边带，同时基频整数值的幅值开始衰减，这说明系统进入强非线性状态。这个判断在齿轮故障诊断里很有价值，因为实际设备中侧隙是无法直接测量的，但可以通过频谱特征反推。

5. 从相图到分岔：参数扫描定位稳定运行区间

5.1 Poincare截面：把连续轨迹离散成周期点

相图看全局，但难以区分周期2和周期4。更精细的作法是取每个激励周期末的状态点，投影到相平面上。比如激励周期 T = 2π/ω，那么在每个 t = kT 时刻记录位移和速度，得到的点集就是Poincare截面。周期1运动在截面上只留下1个点，周期2是2个点，混沌则呈现为分形分布的密集点。这个思想在非线性动力学里是标准工具，拿到MATLAB里实现却很轻量。

% 计算Poincare截面 T_cycle = 2*pi / param.w; k_steps = 100; % 每个周期取100步 n_cycle = 100; % 取100个周期 points = zeros(n_cycle, 2); for i = 1:n_cycle t_span = [i-1, i] * T_cycle; [~, z_temp] = ode45(@(t,z) gearODE(t,z,param), t_span, z0); % 用最后一个时刻的状态作为截面点 points(i, :) = z_temp(end, [1 2]); z0 = z_temp(end, :); % 更新初始条件，连续积分 end plot(points(:,1)*1e6, points(:,2), 'b.');

这段代码按周期拆积分，每个周期结束时的状态就是截面点。注意z0要逐周期更新，否则系统状态不连续。如果某个周期解已经稳定，所有点会重叠在一起；如果是周期2，会看到两个分离的点。实际运行时会发现，用固定步长比ode45默认的变步长更稳定，可以在ode45里加options = odeset('MaxStep', T_cycle/100)。

5.2 用转速扫描画分岔图

分岔图是Poincare截面随着某个参数（通常是转速）连续变化时的汇总图。常见做法是让转速从低到高扫描，每个转速下先去掉前若干个周期的瞬态，再保存后几十个周期的位移值，最后把所有点画成散点图。下面给出一个可直接套用的框架：

sweep_rpm = 800:50:3000; % 转速扫描范围，单位rpm num_cycles_skip = 30; % 跳过瞬态周期 num_cycles_save = 50; % 保存稳态周期 bifur_x = []; for rpm = sweep_rpm param.w = 2*pi*rpm/60 * z_teeth; % z_teeth为齿数 z_current = [1e-5; 0]; T_cycle = 2*pi / param.w; % 第一段：丢弃瞬态 [~, z_skip] = ode45(@(t,z) gearODE(t,z,param), ... linspace(0, num_cycles_skip*T_cycle, num_cycles_skip*200), z_current); z_current = z_skip(end, :).'; % 第二段：记录稳态周期点 [~, z_steady] = ode45(@(t,z) gearODE(t,z,param), ... linspace(0, num_cycles_save*T_cycle, num_cycles_save*200), z_current); % 提取每个周期末的位移 for i = 1:num_cycles_save idx = round(i * 200); % 因为用linspace每周期取200步 bifur_x = [bifur_x; z_steady(idx, 1)]; end end figure('Color','w'); plot(repmat(sweep_rpm, num_cycles_save, 1), bifur_x*1e6, 'b.', 'MarkerSize', 1); xlabel('转速 (rpm)'); ylabel('Poincare位移 (μm)');

这段代码用linspace强制均匀步长，让每个周期的采样点数完全一致，这样取idx时才不会错位。如果不做均匀时间划分，每周期末的时间点不一定落在返回的时间数组上，截取就会引入误差。repmat是为了让每个转速下的50个点并排画在对应的横坐标上。

5.3 工程应用：用分岔图避开混沌区

分岔图有一个直接用途：观察系统在哪个转速区间发生周期倍化。比如你从图中看到在1200~1400rpm区域，Poincare点从1个变成2个，说明周期2开始出现；到1600rpm以后点云弥散，对应混沌振动。设计时就要尽量让工作转速避开这些区域，或者通过加大阻尼比把混沌窗口压缩。

更实用的技巧是结合第4章的FFT，把分岔图和最大幅值曲线对照看。当分岔图中出现混沌时，频谱会出现宽带噪声基底，同时啮合频率处的能量明显向边带转移。这个特征比单纯的时域振幅更能说明问题。日常调试中，我一般先用本章的扫描方法快速定位可疑转速区间，再做一次高频采样FFT确认频率结构，这个方法比盲目跑有限元省太多时间。

上面给出的代码片段里，odeset的MaxStep设置很关键。如果步长太大，间隙非线性处的突变会被平滑，系统可能人为地显示出更强的周期性；步长太小则计算时间成倍增加，一般取激励周期的1/100到1/200之间即可。我在跑不同齿数参数时发现，齿数越多，啮合频率越高，必须相应减小MaxStep才能保持一样的分辨率。验证你的分岔图是否正确有一个土办法：取分岔图上看起来最乱的那个转速，单独画出该转速下的Poincare截面，如果截面上的点呈现出清晰的层叠结构而不是完全随机，说明计算是对的，系统确实处于高维混沌状态。这时候再对比相图，你会发现轨迹始终被约束在一个有限区域内——那正是齿轮非线性动力学里最有意思的地方。

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