MATLAB轴承动力学代码功能说明文章

一、代码整体框架与核心目标

本套MATLAB代码围绕6205深沟球轴承展开，基于滚动轴承故障机理构建动力学模型，通过ODE45求解器完成数值计算，可实现正常、内圈故障、外圈故障、滚动体故障四种工况的仿真。核心目标是输出多维度振动数据与力学特征，为后续故障诊断算法验证、剩余寿命（RUL）预测模型训练提供基础数据支撑。

代码整体分为两类核心文件，一是动力学方程定义文件（如Bearing_Normal.m、Bearing_Fault_In.m等），负责构建不同故障模式下的力学模型；二是求解与结果可视化文件（如Solve_Normal.m、Solve_Fault_Out.m等），负责调用求解器、处理数据并输出图表，同时包含FFT.m工具文件用于频谱分析。

二、故障机理与数学模型构建

2.1 四自由度动力学方程

代码采用四自由度模型（内圈x/y方向、外圈x/y方向），忽略z方向运动（文档4答复10），基本动力学方程形式如下，各参数含义及取值参考西储大学6205轴承标准参数（文档2表格）：

• 内圈x方向：$m_i{\ddot{x}}_i+c_i{\dot{x}}_i+k_i{x}_i=-F_x+m_{i}e{\omega }^2\cos (\omega t)+W_x$
• 内圈y方向：$m_i{\ddot{y}}_i+c_i{\dot{y}}_i+k_i{y}_i=-F_y+m_{i}e{\omega }^2\sin (\omega t)-m_{i}g+W_y$
• 外圈x方向：$m_o{\ddot{x}}_o+c_o{\dot{x}}_o+k_o{x}_o=F_x$
• 外圈y方向：$m_o{\ddot{y}}_o+c_o{\dot{y}}_o+k_o{y}_o=F_y-m_{o}g$

其中，$m_i/m_o$为内外圈质量，$c_i/c_o$为支撑阻尼，$k_i/k_o$为支撑刚度，$F_x/F_y$为滚动体与滚道的总接触力，$e$为内圈偏心量，$\omega$为角速度，$W_x/W_y$为径向载荷，$g$为重力加速度。

2.2 关键力学模型推导

（1）接触刚度计算

基于赫兹接触理论，球轴承接触刚度需通过主曲率参数推导（文档2、3），核心公式如下：

1. 主曲率和计算：

   • 内圈接触主曲率和：${\rho }_{in}=\frac{2}{D_b}+\frac{2}{D_b}+\frac{2}{D_m-D_b}-\frac{1}{f_i{D}_b}$
   • 外圈接触主曲率和：${\rho }_{out}=\frac{2}{D_b}+\frac{2}{D_b}-\frac{2}{D_m+D_b}-\frac{1}{f_o{D}_b}$
     （$D_m$为节圆直径，$D_b$为滚动体直径，$f_i/f_o$为内外圈沟道曲率系数，取值0.515/0.525）

2. 单对接触刚度：
   $k_{bi}=\frac{2\sqrt{2}}{3}\cdot \frac{E}{1-{\mu }^2}\cdot {\delta }_i^{-3/2}\cdot {\rho }_{in}^{-1/2}$
   $k_{bo}=\frac{2\sqrt{2}}{3}\cdot \frac{E}{1-{\mu }^2}\cdot {\delta }_o^{-3/2}\cdot {\rho }_{out}^{-1/2}$
   （$E$为弹性模量，取2.07e11 Pa；$\mu$为泊松比，取0.3；${\delta }_i/{\delta }_o$为接触位移系数，取0.577/0.68）

3. 等效接触刚度：
   $k_b={\left(\frac{1}{(1/k_{bi}{)}^{2/3}+(1/k_{bo}{)}^{2/3}}\right)}^{3/2}$

（2）时变位移激励（故障引入）

不同故障通过“时变位移激励”$h$引入模型，模拟缺陷对接触的冲击（文档2、4）：

• 外圈故障：$h_{out}$由外圈缺陷角度${\varphi }_{os}$、缺陷宽度对应的角度${\varphi }_{do}$控制，仅当滚动体运动至缺陷区域时产生激励。
• 内圈故障：$h_{in}$与内圈转速关联，缺陷角度随内圈转动实时更新（${\varphi }_{is}=\text{mod}(\omega t+{\beta }_0,2\pi )$，${\beta }_0$为初始缺陷角）。
• 滚动体故障：$h_{ball}$由滚动体自旋角度${\varphi }_b$控制，缺陷对称分布（${\varphi }_b=\text{mod}({\omega }_{b}t,2\pi )$，${\omega }_b$为滚动体自转角速度）。

（3）接触力计算

1. 接触变形量：$\delta (j)=(x_i-x_o)\cos \theta (j)+(y_i-y_o)\sin \theta (j)-c_d-h$
   （$\theta (j)$为第j个滚动体的方位角，$c_d$为均化后的径向游隙，$h$为时变位移激励）

2. 亥维塞函数（接触开关）：
   $H(j)=\{\begin{array}{ll}1& \delta (j)>0\text{（接触，产生力）}\\ 0& \delta (j)\le 0\text{（无接触，力为0）}\end{array}$

3. 单滚动体接触力：
   $F_x(j)=k_b\cdot \delta (j{)}^{3/2}\cdot H(j)\cdot \cos \theta (j)$
   $F_y(j)=k_b\cdot \delta (j{)}^{3/2}\cdot H(j)\cdot \sin \theta (j)$

4. 总接触力：$F_x={\sum }_{j=1}^{N_b}{F}_x(j)$，$F_y={\sum }_{j=1}^{N_b}{F}_y(j)$（$N_b$为滚动体个数，取9）

2.3 特征频率验证

代码内置各故障特征频率计算，确保仿真数据与理论一致（文档2、3），以1000 rpm为例：

• 转频：$f_r=n/60=16.67$ Hz
• 外圈故障频率：$f_{BPFO}=\frac{N_{b}n}{2\times 60}\cdot (1-\frac{D_b}{D_m})=59.73$ Hz
• 内圈故障频率：$f_{BPFI}=\frac{N_{b}n}{2\times 60}\cdot (1+\frac{D_b}{D_m})=90.27$ Hz
• 滚动体故障频率：$f_{BSF}=\frac{D_{m}n}{2\times 60\times D_b}\cdot [1-(\frac{D_b}{D_m}{)}^2]=39.24$ Hz

三、核心代码模块功能解析

3.1 动力学方程定义文件（以Bearing_Fault_Out.m为例）

该类文件为ODE45求解器提供微分方程函数，核心流程如下：

1. 参数初始化：定义轴承几何参数（D、d、Dm等）、材料参数（E、μ）、工况参数（n、Wx）、故障参数（缺陷宽度L）。
2. 接触刚度计算：按2.2.1节公式计算$k_{bi}$、$k_{bo}$、$k_b$。
3. 时变激励与接触力计算：循环计算每个滚动体的方位角$\theta (j)$、接触变形量$\delta (j)$，通过亥维塞函数判断接触状态，最终求和得到总接触力$F_x/F_y$。
4. 微分方程输出：将四自由度动力学方程转化为一阶微分方程组（8个状态量：x_i, ${\dot{x}}_i$, y_i, ${\dot{y}}_i$, x_o, ${\dot{x}}_o$, y_o, ${\dot{y}}_o$），输出dy向量。

3.2 求解与可视化文件（以Solve_Fault_In.m为例）

该类文件负责调用求解器并输出结果，核心流程如下：

1. 求解器配置：设置时间步长（tspan1=1e-5）、仿真时长（0-6s）、初始状态（y0，微小位移与零速度）、求解精度（RelTol=1e-3，AbsTol=1e-6）。
2. 数值求解：调用ode45求解微分方程，截取前0.5s过渡段后的数据（避免初始扰动影响）。
3. 数据处理：通过diff函数计算加速度（dy = diff(y)./diff(t)），补充末端数据保证维度一致。
4. 结果可视化：生成7类核心图表，具体如下表：

图表编号	图表类型	核心输出内容	物理意义
1	振动位移时域图	内圈x/y、外圈x/y方向位移（5-5.05s局部放大）	观察故障冲击导致的位移突变
2	振动速度时域图	内圈x/y、外圈x/y方向速度	反映振动能量变化
3	相图	位移-速度曲线（1-3s数据）	判断系统稳定性，故障时曲线更紊乱
4	轴心轨迹	内圈x-y、外圈x-y位移曲线	正常轴承为闭合轨迹，故障时轨迹不规则
5	振动加速度图	内圈x/y、外圈x/y方向加速度（4.75-5.25s放大）	故障冲击在加速度信号中更明显
6	频谱与包络谱	加速度信号的频谱（0-6000Hz）、包络谱（0-500Hz）	提取故障特征频率（如BPFI、边频带）
7	滚道接触力图	1-4号滚动体的接触力随时间变化	观察故障对接触力分布的影响

3.3 工具文件（FFT.m）

实现快速傅里叶变换，将时域信号转化为频域信号，核心步骤：

1. 计算采样频率$fs=\text{size}(t,1)/(t_{\text{end}}-t_{\text{start}})$。
2. 对信号进行FFT变换并归一化，取正频率部分（0-fs/2）。
3. 修正幅值（除直流分量外，其他频率幅值乘以2），输出频率向量f与幅值向量p。

四、参数可调性与扩展应用

4.1 参数调整指南

代码支持多维度参数修改，关键可调参数及调整位置如下：

• 工况参数：转速n（动力学方程文件中）、径向载荷Wx/Wy（动力学方程文件中）、仿真时长tspan（求解文件中）。
• 故障参数：缺陷宽度L（动力学方程文件中）、外圈缺陷方位（外圈故障文件中${\varphi }_{os}$，对应4-8点钟方向）。
• 轴承参数：更换轴承型号时，需修改D、d、Dm、Db、Nb等几何参数（所有文件保持一致），滚子轴承需额外修改接触刚度公式（文档4答复15，接触穿透指数由3/2改为10/9）。

4.2 扩展应用场景

1. 故障诊断算法验证：输出的频谱、包络谱可用于验证故障特征提取算法（如小波变换、经验模态分解EMD），时域数据可用于训练智能诊断模型（如CNN、LSTM）。
2. 剩余寿命（RUL）预测：需补充缺陷扩展模型（文档4提问19），将固定缺陷宽度L改为随时间增长的函数（如$L(t)=L_0+vt$，v为扩展速率），生成全生命周期振动数据，用于RUL预测模型训练。
3. 多数据集适配：修改轴承参数（如节圆直径Dm、滚动体直径Db）后，可适配XJTU-SY、PHM2012等公开数据集的轴承类型（文档4提问8、17）。

五、仿真结果验证与注意事项

5.1 结果验证标准

参考文档2的故障数值计算效果，正常轴承2倍转频误差约4.02%，外圈故障频率误差约0.13%，内圈故障频率误差约0.097%，滚动体故障频率误差约0.084%，仿真结果精度满足工程需求。

5.2 常见问题与解决方法

1. 频谱中找不到转频/故障频率：检查转速n是否与理论计算一致（文档4问题3、4），节圆直径Dm需与轴承型号匹配（滚动体故障需将39.0e-3改为38.05e-3，文档4问题5）。
2. 缺陷宽度修改后报错：缺陷宽度L需小于滚动体直径Db（文档4问题6），建议取值范围0.1-0.5mm（单位为米，代码中需用科学计数法，如0.5e-3）。
3. 滚动体故障频率与理论不符：滚动体可能同时与内外圈接触，导致频率出现2倍关系（文档4问题23），70Hz与140Hz均为合理结果，需结合具体工况判断。