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多行的odefunfunction 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/21); 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运动窗口对应着现场测试中齿轮箱异响的工况后来被证明是间隙参数设置过大的后果。最后展示下这个模型的威力——当输入扭矩突变时庞加莱图上的极限环突然发散对应着轮齿的瞬时脱啮。这种非线性现象用传统的频响分析法根本捕捉不到非得靠时域仿真才能揪出来。