代码实现了一个完整的滚动轴承故障动力学仿真系统，主要通过对轴承系统的多体动力学建模来模拟不同故障状态下的振动响应特性。

系统首先基于轴承的几何参数和材料属性，利用赫兹接触理论计算内外圈与滚动体之间的接触刚度系数，建立轴承系统的非线性动力学微分方程。在动力学仿真过程中，系统考虑了轴承的多种故障类型，包括外滚道故障、内滚道故障和滚动体故障，通过引入相应的故障位移激励函数来模拟实际故障对轴承动态行为的影响。

系统采用数值积分方法求解包含故障激励的非线性微分方程组，获得轴承各部件的位移、速度和加速度响应。通过对仿真得到的振动信号进行时域和频域分析，系统能够揭示不同故障类型在振动信号中表现出的特征模式，为轴承故障诊断提供理论依据和特征提取指导。整个仿真过程涵盖了从轴承参数计算、动力学建模、故障模拟到信号分析的完整流程，为轴承设计和故障诊断研究提供了有效的数值仿真工具。

滚动轴承故障动力学仿真系统流程图
===================================


轴承参数初始化
↓
[几何参数输入] → [材料属性设置] → [工况条件定义]
↓
接触刚度计算
↓
[赫兹接触理论] → [接触变形计算] → [刚度系数拟合]
↓
故障参数设置
↓        ↓        ↓
[外圈故障]  [内圈故障]  [滚动体故障]
↓        ↓        ↓
[故障宽度]  [故障位置]  [故障角度]
↓
动力学方程构建
↓
[多体动力学模型] → [非线性接触力] → [故障激励函数]
↓
数值求解
↓
[微分方程组] → [龙格库塔法] → [时域响应计算]
↓
信号分析
↓        ↓        ↓
[加速度提取]  [时域分析]  [频域分析]
↓
结果输出
↓
[振动信号保存] → [特征图谱生成] → [故障特征识别]

详细算法步骤

第一步进行轴承系统参数初始化，根据轴承型号读取几何尺寸参数，包括滚动体数量、内外圈半径、节圆半径等几何特性，同时设置轴承材料的弹性模量和泊松比等物理参数。

第二步基于赫兹接触理论计算轴承接触刚度，通过理论公式推导内外圈与滚动体之间的接触变形关系，采用曲线拟合方法得到接触力与变形量的非线性刚度系数。

第三步设置故障模拟参数，根据不同类型的故障（外圈故障、内圈故障、滚动体故障）分别定义故障的几何尺寸、位置角度等特征参数，构建相应的故障位移激励函数。

第四步建立轴承系统的多体动力学微分方程，考虑内外圈的质量惯性、轴承支撑刚度阻尼以及滚动体的周期性接触作用，形成包含故障激励的非线性动力学方程组。

第五步采用数值积分方法求解动力学方程，使用变步长龙格库塔法对微分方程组进行数值求解，获得轴承系统在故障状态下的动态响应时间序列。

第六步进行振动信号后处理分析，从求解结果中提取关键测点的加速度响应信号，通过时域波形分析和包络谱分析揭示不同故障类型的特征频率成分。

第七步输出仿真结果和特征图谱，将计算的振动数据保存为文件格式，同时生成直观的时域波形图和频谱图用于故障特征的视觉识别和对比分析。