轴承动力学模型paper复现。 输出质心运动轨迹。 代码注释详细，书写工整规范以便供学者学习轴承动力学。 常量和变量的定义按照论文中的谐音编写，以便学者对照论文和模型进行学习。 因此对小白学习非常友好。 采用自己编写的龙格库塔算法对动力学方程求解。 输出了不同周期下内圈的质心运动轨迹，内圈接触力等动态响应。

轴承动力学模型复现这事儿，说难不难但细节特磨人。咱先看核心方程：内圈质心运动方程长得像俩弹簧拴着铁球，左边是惯性项，右边是接触力加阻尼。重点在于那个非线性赫兹接触力Qc的计算，论文里写成Qc=Kc*δ^1.5，但δ的具体表达式容易搞错符号。

先上硬货——微分方程组定义。注意变量名直接对应论文里的符号：

def dynamic_eq(t, Y):
    x, dx, y, dy = Y
    
    # 赫兹接触变形量（注意绝对值处理）
    delta_x = np.abs(x - Cj*np.sin(Omga_i*t))
    delta_y = np.abs(y - Cj*np.cos(Omga_i*t))
    
    # 非线性接触力计算
    Qcx = Kc * delta_x**1.5 * np.sign(x - Cj*np.sin(Omga_i*t))
    Qcy = Kc * delta_y**1.5 * np.sign(y - Cj*np.cos(Omga_i*t))
    
    # 阻尼力计算（论文式3.7）
    Dx = Cb * dx 
    Dy = Cb * dy
    
    # 方程组拆解（对应论文式3.5）
    d2x = (Qcx - Dx) / Mi
    d2y = (Qcy - Dy - Mi*g)/Mi
    
    return [dx, d2x, dy, d2y]

这里有个坑：计算δ时绝对值和符号函数得分开处理，否则会丢失方向信息。有个博士师兄在这个地方卡了俩礼拜，最后发现是sign函数套在了绝对值外面。

自己撸龙格库塔才够劲，比调库有意思多了。四阶RK的关键在于分阶段计算斜率：

def runge_kutta_4(func, Y0, t):
    n = len(t)
    Y = np.zeros((n, len(Y0)))
    Y[0] = Y0
    
    for i in range(n-1):
        h = t[i+1] - t[i]
        k1 = func(t[i], Y[i])
        k2 = func(t[i] + h/2, Y[i] + h/2*np.array(k1))
        k3 = func(t[i] + h/2, Y[i] + h/2*np.array(k2)) 
        k4 = func(t[i] + h, Y[i] + h*np.array(k3))
        
        Y[i+1] = Y[i] + h/6*(np.array(k1) + 2*np.array(k2) + 2*np.array(k3) + np.array(k4))
    
    return Y

注意这里k1~k4的处理方式，数组运算必须显式转换，否则会引发维度错误。有个细节是时间步长h不固定，所以每次循环都要重新计算，这对变步长算法友好。

轴承动力学模型paper复现。 输出质心运动轨迹。 代码注释详细，书写工整规范以便供学者学习轴承动力学。 常量和变量的定义按照论文中的谐音编写，以便学者对照论文和模型进行学习。 因此对小白学习非常友好。 采用自己编写的龙格库塔算法对动力学方程求解。 输出了不同周期下内圈的质心运动轨迹，内圈接触力等动态响应。

运行参数设置要讲究，特别是无量纲化处理：

# 参数对照论文Table 1
Mi = 0.8    # 内圈质量
Cj = 5e-4   # 间隙量
Kc = 1.3e9  # 赫兹接触刚度 
Cb = 1500   # 阻尼系数
Omga_i = 300 * 2*np.pi / 60  # 转频转换
g = 9.8

转速转换容易翻车，记住1rpm=2π/60 rad/s。有个硕士生把转速直接当角频率用，结果轨迹图变成心电图，当场自闭。

结果可视化时得注意采样间隔，尤其观察周期性：

plt.figure(figsize=(12,6))
plt.plot(t[::10], sol[::10,0], 'b', alpha=0.7, label='X位移')  # 降采样显示
plt.plot(t[::10], sol[::10,2], 'r', alpha=0.7, label='Y位移')
plt.xlabel('时间(s)')
plt.ylabel('位移(m)')
plt.title('内圈质心运动轨迹')
plt.legend()

这里用了::10进行降采样，不然20000个数据点直接绘图会卡成PPT。alpha参数调节透明度能更好显示轨迹重叠部分。

最后说个实战经验：当接触力出现负值时，多半是赫兹接触力的符号处理有误。这时候应该检查delta计算是否带绝对值、sign函数应用是否正确。另外建议用scipy的solve_ivp做结果对比，避免算法实现错误导致整组数据报废。