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/2+1); 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以下），否则混沌区域会失真；庞加莱图的采样时机必须严格对应啮合周期；处理间隙非线性时，初始条件没给对会导致齿轮卡死。后来加了段初始化代码自动计算静平衡位置，这才解决。