MATLAB六自由度齿轮弯扭耦合动力学代码：含时变啮合刚度、齿侧间隙的振动分析与数值计算

MATLAB六自由度齿轮弯扭耦合动力学代码考虑时变啮合刚度、齿侧间隙根据集中质量法建模含数学方程建立和公式推导并在MATLAB中采用ODE45进行数值计算。 输出齿轮水平和竖直方向的振动位移、振动速度、振动加速度、轮齿间动态啮合力、相图、庞加莱图、分岔图、频谱图齿轮系统运转时出现的异常振动常常让工程师头疼。今儿咱们就动手撸一套MATLAB代码把齿轮传动中的弯扭耦合振动给模拟出来。重点在于处理时变啮合刚度和齿侧间隙这两个非线性因素这俩玩意儿直接影响了齿轮的振动特性。先看动力学模型。假设一对齿轮副每个齿轮考虑平移振动x,y方向和扭转振动总共6个自由度。用集中质量法建模时运动方程长得像这样MATLAB六自由度齿轮弯扭耦合动力学代码考虑时变啮合刚度、齿侧间隙根据集中质量法建模含数学方程建立和公式推导并在MATLAB中采用ODE45进行数值计算。 输出齿轮水平和竖直方向的振动位移、振动速度、振动加速度、轮齿间动态啮合力、相图、庞加莱图、分岔图、频谱图Mq Cq K*q F(t)其中刚度矩阵K里藏着时变啮合刚度k(t)这个k(t)可以近似用正弦函数表示k_mesh k0*(1 0.2*sin(omega_mesh*t)); % 时变刚度系数齿侧间隙的处理得用分段函数这里用了个巧妙的平滑处理避免数值震荡function f backlash_func(delta, b) f delta - b*(tanh(1e4*(delta - b)) tanh(1e4*(delta b)))/2; end接着上ODE45求解器。关键是把二阶微分方程转成一阶状态方程。来看代码片段function dydt gear_ode(t,y) % 解包状态变量 x1 y(1); dx1 y(2); y1 y(3); dy1 y(4); theta1 y(5); dtheta1 y(6); % 同理解包齿轮2状态... % 计算相对位移 delta (x2 - x1)*sin(phi) (y2 - y1)*cos(phi) r1*theta1 - r2*theta2; % 考虑齿侧间隙 delta_eff backlash_func(delta, b); % 时变刚度计算 k k0*(1 0.2*sin(omega_mesh*t)); % 非线性啮合力 F_mesh k*delta_eff c_mesh*(dtheta1*r1 - dtheta2*r2); % 组装加速度项 ddx1 (-cx*dx1 - kx*x1 F_mesh*sin(phi))/m1; % 其他加速度项同理... dydt [dx1; ddx1; dy1; ddy1; dtheta1; ddtheta1; ...]; % 共12个状态量 end跑完仿真后振动响应分析才是重头戏。频谱图用FFT搞定Y fft(acceleration); P2 abs(Y/L); P1 P2(1:L/21); P1(2:end-1) 2*P1(2:end-1); f Fs*(0:(L/2))/L; plot(f,P1)分岔图需要参数扫描这里用了个粗暴但有效的循环for omega 1000:50:3000 % 转速扫描范围 [t,y] ode45(gear_ode, [0 2], init_cond); % 取稳态数据 plot(omega, y(end-1000:end,1), .k); hold on end相图和庞加莱图能直观显示系统状态。注意庞加莱截面的取法要和激励频率同步poincare_times 0:2*pi/omega_mesh:max(t); poincare_data interp1(t, y(:,1), poincare_times); plot(poincare_data(100:end), diff(poincare_data(99:end)), .)这套代码跑下来你会明显看到随着转速变化系统从周期运动进入混沌状态——分岔图上那些突然散开的点就是证据。时变刚度导致的参数激励效应在频谱图上表现为边频带而齿侧间隙会引发次谐波共振。调试时遇到过几个坑ODE45的相对误差容限得调小1e-6以下否则混沌区域会失真庞加莱图的采样时机必须严格对应啮合周期处理间隙非线性时初始条件没给对会导致齿轮卡死。后来加了段初始化代码自动计算静平衡位置这才解决。