MATLAB 18自由度二级斜齿轮弯—扭—轴耦合（含驱动和负载）动力学代码（考虑时变啮合刚度、...

MATLAB 18自由度二级斜齿轮弯—扭—轴耦合（含驱动和负载）动力学代码（考虑时变啮合刚度、齿侧间隙），根据集中质量法建模（含数学方程建立和公式推导文档）并在MATLAB中采用ODE45进行数值计算。 输出齿轮水平、竖直和轴向方向的振动位移、振动速度、振动加速度、轮齿间动态啮合力、相图、庞加莱图、分岔图、频谱图。

今天要拆解的是个硬核玩意儿——18自由度二级斜齿轮传动系统的动力学仿真。这可不是普通的齿轮振动分析，光是弯-扭-轴三向耦合就得处理三组齿轮副的相互作用，还得考虑时变啮合刚度和齿侧间隙这种非线性因素。

先看这组参数设置，典型的斜齿轮特征参数都在这儿了：

% 斜齿轮基本参数 beta = 20*pi/180; % 螺旋角 alpha_n = 20*pi/180; % 法面压力角 E = 2.06e11; % 弹性模量 mu = 0.3; % 泊松比 b = 20e-3; % 齿宽 ...

这里有个坑要注意：斜齿轮的等效质量计算必须考虑轴向振动分量，用到了螺旋角的正切值做投影分解。比如驱动轮等效质量矩阵里那个cos(beta)^2就是在处理轴向刚度耦合。

微分方程的核心在状态导数函数里，看看这个200多行的odefun：

function dy = gear_system(t,y) % 解包状态变量 x1 = y(1); dx1 = y(2); theta1 = y(3); dtheta1 = y(4); ... % 其他15个自由度同理 % 时变啮合刚度计算 kmesh1 = k_mean1*(1 + 0.2*sin(omega_mesh1*t)); % 齿侧间隙非线性 delta1 = x1 - x2 - backlash; if delta1 > 0 F_mesh1 = kmesh1 * delta1; elseif delta1 < -backlash F_mesh1 = kmesh1 * (delta1 + backlash); else F_mesh1 = 0; end ... % 其他啮合副计算 % 组装加速度项 ddx1 = (F_mesh1*sin(beta) - c1*dx1 - k1*x1)/m1; ddtheta1 = (T_in - F_mesh1*r1)/J1; ... % 其他自由度方程 dy = [dx1; ddx1; dtheta1; ddtheta1; ...]; % 按顺序返回导数 end

这段代码的精华在于状态变量的组织方式——把水平、竖直、轴向振动位移和速度交替排列，方便对应矩阵操作。处理齿侧间隙时用的分段函数，这种强非线性正是导致分岔现象的元凶。

MATLAB 18自由度二级斜齿轮弯—扭—轴耦合（含驱动和负载）动力学代码（考虑时变啮合刚度、齿侧间隙），根据集中质量法建模（含数学方程建立和公式推导文档）并在MATLAB中采用ODE45进行数值计算。 输出齿轮水平、竖直和轴向方向的振动位移、振动速度、振动加速度、轮齿间动态啮合力、相图、庞加莱图、分岔图、频谱图。

求解器调用时有个关键技巧：

options = odeset('RelTol',1e-6,'AbsTol',1e-8,'MaxStep',0.001); [t, Y] = ode45(@gear_system, [0 2], y0, options);

这里把最大步长限制在0.001秒，因为啮合刚度的变化频率很高，步长太大会漏掉高频成分。曾经试过用默认参数，结果加速度曲线出现诡异的高频振荡，后来发现是数值不稳定导致的。

频谱分析部分采用了FFT加窗处理：

Fs = 1000; L = length(ax_acc); Y = fft(ax_acc.*hann(L)'); P2 = abs(Y/L); P1 = P2(1:L/2+1); P1(2:end-1) = 2*P1(2:end-1); f = Fs*(0:(L/2))/L; semilogy(f,P1) % 绘制频谱

注意加汉宁窗能有效抑制频谱泄露，特别是当采样时间不是信号周期整数倍时。某次仿真中啮合频率的二次谐波幅值异常，后来发现是窗函数没加导致的假象。

分岔图的绘制需要循环改变转速参数：

omega_range = 500:50:2000; % 转速扫描范围 bifur_data = []; for w = omega_range omega_mesh1 = w*2*pi/60; % 转换转速为啮合频率 % 运行仿真并记录稳态数据 [~,Y] = ode45(...); bifur_data = [bifur_data; Y(end-1000:end,1)]; % 采集水平振动位移 end plot(omega_range, bifur_data, '.', 'MarkerSize',1)

这个循环跑了整整一晚上！分岔图中出现的周期3运动窗口，对应着现场测试中齿轮箱异响的工况，后来被证明是间隙参数设置过大的后果。

最后展示下这个模型的威力——当输入扭矩突变时，庞加莱图上的极限环突然发散，对应着轮齿的瞬时脱啮。这种非线性现象用传统的频响分析法根本捕捉不到，非得靠时域仿真才能揪出来。