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运动窗口，对应着现场测试中齿轮箱异响的工况，后来被证明是间隙参数设置过大的后果。

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