MATLAB轴承—转子系统动力学代码（可考虑轴承故障、转子碰摩、转子涂层），根据集中质量法建模（含数学方程建立和公式推导）并在MATLAB中采用ODE45进行数值计算。 可模拟不同系统参数的动力学特性，输出时域加速度、速度和位移波形、动态碰摩力、频谱图、相图、庞加莱图和分岔图。
最近在折腾旋转机械的故障仿真，发现用集中质量法建模轴承-转子系统还挺有意思。这个模型能模拟轴承故障、转子碰摩这些典型问题，今天就来说说怎么用MATLAB的ODE45实现这个动力学仿真。
先上核心的数学模型。我们把转子简化为两个垂直方向的自由度，考虑质量块的位移x和y。轴承刚度用k表示，阻尼系数c，碰摩发生时会产生非线性接触力。系统方程长这样：
mx'' + cx' + k*x = Funbalance + Frubbing
my'' + cy' + k*y = Funbalance + Frubbing
这里碰摩力F_rubbing的模型是关键。当转子位移超过间隙δ时触发碰摩，可以用库伦摩擦模型：
if sqrt(x^2 + y^2) > δ
F_rubbing = -μN(x/r, y/r)
else
F_rubbing = 0
end
接下来是MATLAB的实现部分。先定义系统参数：
m = 12.5; % 质量 kg k = 8e5; % 刚度 N/m c = 1200; % 阻尼 N·s/m mu = 0.08; % 摩擦系数 delta = 1e-4; % 碰摩间隙 m omega = 1200; % 转速 rpm微分方程的函数实现需要注意状态变量的拆分。这里把x和y方向的位移、速度都放进状态向量：
function dydt = rotor_sys(t, y) % y = [x; dx/dt; y; dy/dt] global m c k mu delta omega F_unbalance = 0.05*m*omega^2*sin(omega*t); % 不平衡激励 % 碰摩力计算 r = sqrt(y(1)^2 + y(3)^2); if r > delta F_rub = -mu*k*(r - delta)*[y(1); y(3)]/r; else F_rub = [0; 0]; end % 微分方程组 dydt = zeros(4,1); dydt(1) = y(2); dydt(2) = (-c*y(2) - k*y(1) + F_unbalance + F_rub(1))/m; dydt(3) = y(4); dydt(4) = (-c*y(4) - k*y(3) + F_unbalance + F_rub(2))/m; end这个函数有几个需要注意的点：1）全局变量的使用虽然方便，但后续可以考虑参数传递优化；2）碰摩力的方向处理需要归一化坐标；3）不平衡力的幅值需要根据实际工况调整。
求解器调用采用自适应步长的ODE45：
[t, Y] = ode45(@rotor_sys, [0 2], [1e-5 0 1e-5 0], ... odeset('RelTol',1e-6));结果后处理可以玩出很多花样。时域波形直接plot就行，频谱分析用FFT要注意窗函数选择：
% 加速度计算 acc = gradient(Y(:,2), t(2)-t(1)); % 频谱分析 Fs = 1/(t(2)-t(1)); NFFT = 2^nextpow2(length(acc)); freq = Fs/2*linspace(0,1,NFFT/2); fft_acc = fft(acc, NFFT);庞加莱图通过周期采样实现，分岔图则需要扫频计算。比如改变转速omega，观察系统行为变化：
omega_range = 800:50:2000; bifurcation_data = zeros(length(omega_range), 100); for i = 1:length(omega_range) omega = omega_range(i); [~, Y] = ode45(...); % 重新求解 bifurcation_data(i,:) = Y(end-99:end,1); % 取稳态解 end实际跑代码时可能会发现，当转速接近临界转速时，系统会出现明显的共振峰。加入碰摩后，频谱中会出现超谐波成分，庞加莱图呈现点阵分布，说明系统进入准周期运动状态。
调试过程中遇到过几个坑：1）刚度系数太大导致数值发散，需要调整求解器相对误差；2）碰摩力计算时忘记归一化导致方向错误；3）分岔扫频耗时太长，后来改用parfor并行计算才解决。
最后说个有意思的现象：当摩擦系数超过0.1时，系统会出现混沌运动，相图轨迹不再闭合，分岔图上出现典型的倍周期分岔序列。这说明即使简单的两自由度模型，也能呈现丰富的非线性动力学行为。
