MATLAB实现滚动轴承二自由度动力学建模与故障诊断

1. 滚动轴承动力学分析概述滚动轴承作为旋转机械的核心部件其动力学特性直接影响设备运行稳定性与寿命。在MATLAB中建立二自由度轴承动力学模型能够有效模拟正常状态及各类故障下的动态响应特征。这种仿真方法为轴承状态监测和故障诊断提供了重要技术手段。轴承动力学建模需要考虑以下几个关键因素滚动体与内外圈的弹性接触变形轴承游隙和预紧力润滑条件与摩擦效应局部缺陷引起的激励变化典型的滚动轴承故障包括内圈故障表现为与轴旋转频率相关的周期性冲击外圈故障特征频率与轴承几何参数相关滚动体故障产生与保持架转速相关的调制信号2. 二自由度动力学模型构建2.1 基本运动方程建立采用集中参数法建立轴承的二自由度动力学模型考虑径向(x,y方向)两个自由度mẍ cẋ kx Fx(t) Fd_x(t) mÿ cẏ ky Fy(t) Fd_y(t)其中m等效质量包含轴和轴承质量c阻尼系数考虑油膜阻尼和结构阻尼k时变刚度反映滚动体位置变化F(t)外部载荷Fd(t)缺陷引起的动态激励2.2 轴承刚度计算轴承刚度具有明显的时变特性需考虑滚动体位置变化function K bearing_stiffness(theta, n_balls, k_contact) % theta: 滚动体位置角度(rad) % n_balls: 滚动体数量 % k_contact: 接触刚度 K zeros(2,2); for i 1:n_balls phi_i 2*pi*(i-1)/n_balls; delta_theta theta - phi_i; if abs(delta_theta) contact_angle_threshold K K k_contact * [cos(phi_i)^2, cos(phi_i)*sin(phi_i); cos(phi_i)*sin(phi_i), sin(phi_i)^2]; end end end2.3 故障激励建模不同故障类型对应不同的激励函数内圈故障激励function F inner_race_defect(t, omega, n_balls, D_pitch, d_ball, beta) % omega: 轴旋转速度(rad/s) % D_pitch: 节圆直径 % d_ball: 滚动体直径 % beta: 接触角 f_cage omega/2 * (1 - d_ball/D_pitch * cos(beta)); defect_freq n_balls * f_cage; F defect_amplitude * sin(2*pi*defect_freq*t) .* rectpuls(mod(omega*t,1/defect_freq)); end外圈故障激励function F outer_race_defect(t, omega, n_balls, D_pitch, d_ball, beta) f_cage omega/2 * (1 d_ball/D_pitch * cos(beta)); defect_freq n_balls * f_cage; F defect_amplitude * sin(2*pi*defect_freq*t) .* rectpuls(mod(omega*t,1/defect_freq)); end3. MATLAB求解实现3.1 ODE求解器配置采用ode45求解器处理非线性时变系统% 参数设置 params.m 5; % 等效质量(kg) params.c 500; % 阻尼系数(N·s/m) params.n_balls 8; % 滚动体数量 params.k_contact 1e8; % 接触刚度(N/m) % 初始条件 x0 [0; 0; 0; 0]; % [x; y; dxdt; dydt] % 时间范围 tspan [0 0.1]; % 求解选项设置 options odeset(RelTol,1e-6,AbsTol,1e-8); % 调用求解器 [t, X] ode45((t,x) bearing_equations(t,x,params), tspan, x0, options);3.2 状态方程函数function dxdt bearing_equations(t, x, params) % 状态变量分解 pos x(1:2); vel x(3:4); % 轴承当前角度 theta mod(params.omega * t, 2*pi); % 计算时变刚度 K bearing_stiffness(theta, params.n_balls, params.k_contact); % 故障激励 F_defect [0; 0]; if params.has_inner_defect F_defect F_defect inner_race_defect(t, params); end if params.has_outer_defect F_defect F_defect outer_race_defect(t, params); end % 外部载荷 F_ext [params.Fx; params.Fy]; % 加速度计算 accel params.M \ (F_ext F_defect - params.C*vel - K*pos); % 状态导数 dxdt [vel; accel]; end4. 故障特征分析与可视化4.1 时域响应分析figure; subplot(2,1,1); plot(t, X(:,1), b, LineWidth, 1.5); xlabel(Time (s)); ylabel(Displacement x (m)); title(X方向位移时程); subplot(2,1,2); plot(t, X(:,2), r, LineWidth, 1.5); xlabel(Time (s)); ylabel(Displacement y (m)); title(Y方向位移时程);4.2 频域特征提取% 计算FFT Fs 1/(t(2)-t(1)); % 采样频率 L length(t); % 信号长度 Yx fft(X(:,1)); P2x abs(Yx/L); P1x P2x(1:L/21); P1x(2:end-1) 2*P1x(2:end-1); Yy fft(X(:,2)); P2y abs(Yy/L); P1y P2y(1:L/21); P1y(2:end-1) 2*P1y(2:end-1); f Fs*(0:(L/2))/L; figure; subplot(2,1,1); plot(f, P1x, b, LineWidth, 1.5); xlabel(Frequency (Hz)); ylabel(Amplitude); title(X方向频谱); subplot(2,1,2); plot(f, P1y, r, LineWidth, 1.5); xlabel(Frequency (Hz)); ylabel(Amplitude); title(Y方向频谱);4.3 包络分析用于故障诊断% 希尔伯特变换提取包络 analytic_signal hilbert(X(:,1)); envelope abs(analytic_signal); % 包络谱分析 Y_env fft(envelope); P2_env abs(Y_env/L); P1_env P2_env(1:L/21); P1_env(2:end-1) 2*P1_env(2:end-1); figure; plot(f, P1_env, m, LineWidth, 1.5); xlabel(Frequency (Hz)); ylabel(Amplitude); title(包络谱);5. 不同故障状态的动态响应对比5.1 正常状态特征正常轴承的振动信号具有以下特点时域波形呈现准周期性频谱中以轴频及其谐波为主高频区域能量较低包络谱无明显特征频率5.2 内圈故障特征内圈故障的典型表现时域出现周期性冲击频谱中出现内圈故障特征频率及其边带包络谱中清晰显示故障特征频率冲击间隔与轴转速相关5.3 外圈故障特征外圈故障的识别特征时域冲击间隔固定频谱中外圈故障频率成分突出包络谱中外圈故障频率明显可能伴随高阶谐波5.4 滚动体故障特征滚动体故障的振动特性冲击间隔与保持架转速相关频谱中出现滚动体故障频率包络谱中滚动体通过频率明显常伴有调制边带6. 模型验证与参数影响分析6.1 模型验证方法理论频率验证比较仿真得到的故障特征频率与理论计算值能量分布验证检查各频段能量分布是否符合预期冲击特性验证分析时域冲击间隔与理论值的吻合度% 计算理论故障频率 BPFI n_balls/2 * omega/(2*pi) * (1 d_ball/D_pitch * cos(beta)); % 内圈故障 BPFO n_balls/2 * omega/(2*pi) * (1 - d_ball/D_pitch * cos(beta)); % 外圈故障 FTF omega/(2*pi) * 0.5 * (1 - d_ball/D_pitch * cos(beta)); % 保持架频率 BSF D_pitch/d_ball * omega/(2*pi) * 0.5 * (1 - (d_ball/D_pitch * cos(beta))^2); % 滚动体故障6.2 关键参数影响分析接触刚度影响刚度增大→振动幅值减小高频成分增加阻尼系数影响阻尼增大→振动衰减加快共振峰幅值降低缺陷尺寸影响缺陷越大→冲击幅值越大谐波成分更丰富转速影响转速提高→故障频率线性增加冲击间隔缩短7. 工程应用与扩展7.1 故障诊断系统集成将轴承模型集成到状态监测系统中function [fault_type, severity] diagnose_bearing(signal, fs) % 特征提取 features extract_features(signal, fs); % 加载预训练的分类模型 load(bearing_classifier.mat); % 故障分类 fault_type predict(classifier, features); % 严重程度评估 severity assess_severity(features); end7.2 模型扩展方向考虑非线性因素赫兹接触非线性间隙非线性润滑状态变化多故障耦合分析同时存在多种故障的相互作用故障之间的调制效应不确定性分析参数不确定性影响随机激励的影响智能诊断算法深度学习特征提取迁移学习应用数字孪生技术提示在实际应用中建议先通过理论计算确定轴承各故障特征频率再结合仿真结果进行对比分析。同时要注意实测信号通常包含噪声干扰需要适当增加信号处理环节。