MATLAB轴承动力学代码（正常、外圈故障、内圈故障、滚动体故障），根据滚动轴承故障机理建模（含数学方程建立和公式推导）并在MATLAB中采用ODE45进行数值计算。 可模拟不同轴承故障类型，输出时域加速度波形、滚道接触力、相图。

在机械工程领域，滚动轴承作为关键部件，其运行状态的监测与故障诊断至关重要。今天咱们就聊聊如何利用 MATLAB 对轴承动力学进行模拟，涵盖正常状态以及外圈故障、内圈故障和滚动体故障这几种常见情况。

一、滚动轴承故障机理建模

1. 数学方程建立

滚动轴承动力学模型通常基于牛顿第二定律。以一个简单的单自由度模型为例，设轴承系统的质量为 \(m\)，刚度为 \(k\)，阻尼为 \(c\)，作用在轴承上的外部激励力为 \(F(t)\)。那么，系统的运动方程可以表示为：

\[ m\ddot{x} + c\dot{x} + kx = F(t) \]

其中 \(x\) 是位移，\(\dot{x}\) 是速度，\(\ddot{x}\) 是加速度。

对于滚动轴承不同部位的故障，需要在这个基本方程上进行拓展。比如，当存在外圈故障时，故障点与滚动体的接触会产生周期性的冲击激励，这个激励可以用一个与故障特征频率相关的函数来表示。假设外圈故障特征频率为 \(f{o}\)，则激励力 \(F{o}(t)\) 可以写成：

\[ F{o}(t) = A{o}\sum{n = -\infty}^{\infty}\delta(t - nT{o}) \]

这里 \(A{o}\) 是激励幅值，\(\delta(t)\) 是狄拉克函数，\(T{o} = 1 / f_{o}\) 是故障周期。

同理，内圈故障和滚动体故障也可以类似地通过建立与各自故障特征频率相关的激励力函数来完善运动方程。

2. 公式推导

以求解位移 \(x(t)\) 为例，对于上述运动方程，我们可以将其转化为一阶微分方程组。令 \(y{1} = x\)，\(y{2} = \dot{x}\)，则原方程可写成：

\[ \begin{cases} \dot{y}{1} = y{2} \\ \dot{y}{2} = \frac{1}{m}(F(t) - c y{2} - k y_{1}) \end{cases} \]

这样就把一个二阶微分方程转化为了一阶微分方程组，便于在 MATLAB 中使用数值方法求解。

二、MATLAB 中的数值计算

在 MATLAB 里，我们使用 ode45 函数来进行数值计算。ode45 是基于龙格 - 库塔（Runge - Kutta）算法的一种数值求解器，适用于大多数非刚性问题。

下面是一段简单的 MATLAB 代码示例，用于求解上述单自由度系统的运动方程（假设为正常状态，激励力 \(F(t)\) 为常数 \(F_{0}\)）：

% 参数设置
m = 1; % 质量
c = 0.1; % 阻尼
k = 10; % 刚度
F0 = 1; % 激励力

% 定义时间范围
tspan = 0:0.01:10; 

% 初始条件
y0 = [0; 0]; 

% 定义微分方程
odefun = @(t,y) [y(2); (F0 - c*y(2) - k*y(1))/m]; 

% 使用 ode45 求解
[t,y] = ode45(odefun, tspan, y0); 

% 提取位移和加速度
x = y(:,1);
acceleration = (F0 - c*y(:,2) - k*y(:,1))/m; 

% 绘制位移随时间变化曲线
figure;
plot(t, x);
xlabel('时间 (s)');
ylabel('位移 (m)');
title('正常状态下轴承位移随时间变化');

在这段代码中，首先定义了系统的参数，包括质量、阻尼、刚度和激励力。然后设定了时间范围 tspan 和初始条件 y0。接下来，使用匿名函数 odefun 定义了一阶微分方程组。最后，通过 ode45 函数求解该方程组，并提取位移和加速度信息，绘制出位移随时间变化的曲线。

当处理不同故障类型时，只需要修改激励力函数 odefun 中 F(t) 的定义即可。比如对于外圈故障，可按如下方式修改：

% 参数设置（假设外圈故障相关参数）
m = 1; 
c = 0.1; 
k = 10; 
A_o = 1; % 外圈故障激励幅值
f_o = 10; % 外圈故障特征频率
T_o = 1/f_o; 

% 定义时间范围
tspan = 0:0.01:10; 

% 初始条件
y0 = [0; 0]; 

% 定义微分方程（考虑外圈故障激励）
odefun = @(t,y) [y(2); (A_o*sum(dirac(t - 0:1/f_o:10)) - c*y(2) - k*y(1))/m]; 

% 使用 ode45 求解
[t,y] = ode45(odefun, tspan, y0); 

% 提取位移和加速度
x = y(:,1);
acceleration = (A_o*sum(dirac(t - 0:1/f_o:10)) - c*y(2) - k*y(1))/m; 

% 绘制时域加速度波形
figure;
plot(t, acceleration);
xlabel('时间 (s)');
ylabel('加速度 (m/s^2)');
title('外圈故障时轴承时域加速度波形');

这里通过修改激励力部分，考虑了外圈故障的影响，并绘制出了时域加速度波形。

三、结果输出

通过上述代码，我们可以模拟不同轴承故障类型，并输出时域加速度波形、滚道接触力、相图等关键信息。

1. 时域加速度波形

如前面代码中所示，通过绘制加速度随时间的变化曲线，我们可以直观地观察到不同故障类型下加速度的变化特征。正常状态下，加速度可能相对平稳，而故障状态下，会出现明显的冲击特征。

2. 滚道接触力

滚道接触力可以通过对运动方程中与接触相关的力项进行分析和计算得到。假设接触力为 \(F_{contact}\)，它与位移、速度等状态变量有关，可以通过类似下面的方式计算（这里只是简单示意，实际计算可能更复杂）：

\[ F_{contact} = kx + c\dot{x} \]

在 MATLAB 代码中，可以在求解完微分方程后，根据得到的位移和速度数据进行计算并绘制接触力随时间变化的曲线。

% 计算滚道接触力
F_contact = k*y(:,1) + c*y(:,2); 

% 绘制滚道接触力随时间变化曲线
figure;
plot(t, F_contact);
xlabel('时间 (s)');
ylabel('滚道接触力 (N)');
title('滚道接触力随时间变化');

3. 相图

相图是指速度与位移的关系图，可以帮助我们更深入地理解系统的动力学特性。在 MATLAB 中绘制相图非常简单，只需将位移和速度数据绘制成二维散点图即可。

% 绘制相图
figure;
plot(y(:,1), y(:,2));
xlabel('位移 (m)');
ylabel('速度 (m/s)');
title('轴承运动相图');

通过以上步骤，我们利用 MATLAB 实现了滚动轴承动力学的模拟，从故障机理建模到数值计算，再到关键结果的输出，为滚动轴承的故障诊断提供了有效的分析手段。希望这篇博文对正在研究相关领域的朋友们有所帮助。