1. 轴承振动信号模拟：从理论到实践的完整实现
轴承作为旋转机械的核心部件，其健康状态直接影响设备运行安全。传统故障诊断依赖现场采集数据，但实际工况中故障样本稀缺且获取成本高。通过MATLAB构建轴承振动仿真系统，我们可以在实验室环境下生成各类故障特征的振动信号，为算法验证提供可控数据源。
我曾在某风电齿轮箱诊断项目中，因现场故障样本不足导致算法验证受阻，最终通过仿真系统生成200组不同损伤程度的轴承信号，顺利完成算法测试。这套方法的核心在于将机械动力学原理转化为可计算的数学模型，再通过数值方法求解。下面将详细解析从动力学方程构建到故障特征注入的全流程实现。
2. 动力学模型构建基础
2.1 轴承-转子系统力学建模
健康的滚动轴承振动主要来自滚动体周期性接触产生的激励，其动力学方程可表示为：
M*d2x/dt2 + C*dx/dt + K*x = F(t)其中：
- M为系统质量矩阵（包含内圈、外圈、滚动体质量）
- C为阻尼矩阵（通常取临界阻尼的2%-5%）
- K为刚度矩阵（考虑Hertz接触刚度）
- F(t)为时变激励力
对于6205深沟球轴承的典型参数示例：
m_ball = 0.01; % 单球质量(kg) k_contact = 7e8; % 接触刚度(N/m) c_damping = 1500; % 阻尼系数(N·s/m)2.2 故障特征频率计算
轴承各部件故障特征频率计算公式如下：
| 故障类型 | 计算公式 | 参数说明 |
|---|---|---|
| 内圈故障 | BPFI = (n/2)*fr*(1+d/D*cosα) | n:球数, fr:转频, d:球径, D:节径 |
| 外圈故障 | BPFO = (n/2)*fr*(1-d/D*cosα) | α:接触角 |
| 滚动体故障 | BSF = (D/d)*fr*(1-(d/D*cosα)^2) | |
| 保持架故障 | FTF = (fr/2)*(1-d/D*cosα) |
MATLAB实现示例：
function [BPFI, BPFO, BSF, FTF] = bearing_fault_freq(n, fr, d, D, alpha) BPFI = (n/2)*fr*(1 + d/D*cosd(alpha)); BPFO = (n/2)*fr*(1 - d/D*cosd(alpha)); BSF = (D/d)*fr*(1 - (d/D*cosd(alpha))^2); FTF = (fr/2)*(1 - d/D*cosd(alpha)); end3. 故障激励建模与数值求解
3.1 典型故障激励函数设计
外圈局部损伤会产生周期性冲击，其激励力可建模为：
function F = outer_race_fault(t, BPFO, A) tau = 1/BPFO; % 冲击间隔 Td = 0.001; % 冲击持续时间 F = A * exp(-50*(mod(t,tau)/Td).^2) .* sin(2*pi*3000*t); end内圈损伤因存在载荷周期调制效应，需增加幅值调制项：
function F = inner_race_fault(t, BPFI, A, fr) tau = 1/BPFI; Td = 0.0008; carrier = A * exp(-60*(mod(t,tau)/Td).^2); modulator = 1 + 0.5*sin(2*pi*fr*t); % 转频调制 F = carrier .* modulator .* sin(2*pi*5000*t); end3.2 龙格-库塔数值解法实现
采用ode45求解器进行数值积分时，需将二阶微分方程转化为一阶方程组：
function dxdt = bearing_ode(t, x, M, C, K, F_fn) % 状态变量: x(1)=位移, x(2)=速度 dxdt = zeros(2,1); dxdt(1) = x(2); dxdt(2) = M \ (F_fn(t) - C*x(2) - K*x(1)); end % 调用示例 [t, X] = ode45(@(t,x) bearing_ode(t,x,M,C,K,@(t)outer_race_fault(t,BPFO,10)), ... [0 0.1], [0 0], odeset('RelTol',1e-6));关键参数选择经验：相对误差RelTol建议设为1e-6，时间步长取最高故障频率的1/20（如10kHz采样对应5e-5s步长）
4. 仿真信号后处理与特征增强
4.1 噪声注入与带通滤波
实际信号必然包含噪声，需添加高斯白噪声和工频干扰：
fs = 20e3; % 采样率20kHz noisy_signal = x + 0.02*randn(size(x)) + 0.1*sin(2*pi*50*t'); % 设计80-3000Hz带通滤波器 [b,a] = butter(4, [80 3000]/(fs/2), 'bandpass'); filtered_signal = filtfilt(b, a, noisy_signal);4.2 包络解调分析实现
包络分析是轴承故障诊断的核心方法，MATLAB实现步骤：
- Hilbert变换提取解析信号
- 取模得到包络线
- 对包络谱进行FFT分析
analytic_signal = hilbert(filtered_signal); envelope = abs(analytic_signal); [P_env, f_env] = pwelch(envelope, hamming(1024), 512, 4096, fs); figure; subplot(211), plot(t(1:1000), filtered_signal(1:1000)); subplot(212), plot(f_env, 10*log10(P_env)); xline([BPFO BPFI BSF], 'r--'); % 标注特征频率5. 工程实践中的关键问题
5.1 参数敏感性分析
通过蒙特卡洛仿真评估参数影响（示例仅展示刚度影响）：
k_range = linspace(5e8, 9e8, 20); fault_level = zeros(size(k_range)); for i = 1:length(k_range) [t,X] = ode45(@(t,x)bearing_ode(t,x,M,C,k_range(i),@(t)inner_race_fault(t,BPFI,8,fr)),...); env = abs(hilbert(X(:,1))); P = pwelch(env, [], [], [], fs); fault_level(i) = max(P(f_env>BPFI*0.9 & f_env<BPFI*1.1)); end plot(k_range, fault_level); xlabel('接触刚度(N/m)'); ylabel('故障特征幅值');5.2 常见问题排查指南
| 现象 | 可能原因 | 解决方案 |
|---|---|---|
| 解调谱峰值模糊 | 阻尼系数过大 | 调整阻尼比为临界阻尼的1%-3% |
| 特征频率偏移 | 滚动体打滑率未考虑 | 在BPFI/BPFO公式中添加5%-8%滑差 |
| 高频共振不明显 | 激励力中心频率过低 | 将冲击正弦分量提高到5kHz以上 |
| 数值解发散 | 步长过大或刚度矩阵病态 | 使用ode15s stiff求解器 |
6. 仿真系统进阶扩展
6.1 多故障耦合建模
实际轴承常存在复合故障，需建立耦合激励模型：
function F = compound_fault(t, BPFI, BPFO, A) tau_inner = 1/BPFI; tau_outer = 1/BPFO; F_inner = 0.7*A * exp(-50*(mod(t,tau_inner)/0.001).^2) .* sin(2*pi*5000*t); F_outer = 0.3*A * exp(-60*(mod(t,tau_outer)/0.0012).^2) .* sin(2*pi*4500*t); F = F_inner + F_outer; end6.2 基于Simulink的实时仿真
对于复杂系统，可构建Simulink模型实现硬件在环测试：
- 使用S-Function嵌入动力学方程
- 通过Signal Builder模块注入故障
- 利用DSP System Toolbox进行实时滤波
- 通过xPC Target实现实时仿真
模型关键配置：
- 固定步长ode4（Runge-Kutta）
- 采样时间≤1e-5s
- 启用零交叉检测
