轴承动力学模型paper复现。 输出质心运动轨迹。 代码注释详细书写工整规范以便供学者学习轴承动力学。 常量和变量的定义按照论文中的谐音编写以便学者对照论文和模型进行学习。 因此对小白学习非常友好。 采用自己编写的龙格库塔算法对动力学方程求解。 输出了不同周期下内圈的质心运动轨迹内圈接触力等动态响应。轴承动力学模型复现这事儿说难不难但细节特磨人。咱先看核心方程内圈质心运动方程长得像俩弹簧拴着铁球左边是惯性项右边是接触力加阻尼。重点在于那个非线性赫兹接触力Qc的计算论文里写成QcKc*δ^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[i1] - 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[i1] 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转速转换容易翻车记住1rpm2π/60 rad/s。有个硕士生把转速直接当角频率用结果轨迹图变成心电图当场自闭。结果可视化时得注意采样间隔尤其观察周期性plt.figure(figsize(12,6)) plt.plot(t[::10], sol[::10,0], b, alpha0.7, labelX位移) # 降采样显示 plt.plot(t[::10], sol[::10,2], r, alpha0.7, labelY位移) plt.xlabel(时间(s)) plt.ylabel(位移(m)) plt.title(内圈质心运动轨迹) plt.legend()这里用了::10进行降采样不然20000个数据点直接绘图会卡成PPT。alpha参数调节透明度能更好显示轨迹重叠部分。最后说个实战经验当接触力出现负值时多半是赫兹接触力的符号处理有误。这时候应该检查delta计算是否带绝对值、sign函数应用是否正确。另外建议用scipy的solve_ivp做结果对比避免算法实现错误导致整组数据报废。
基于轴承动力学模型的paper复现：质心运动轨迹输出及龙格库塔算法的编码实践
轴承动力学模型paper复现。 输出质心运动轨迹。 代码注释详细书写工整规范以便供学者学习轴承动力学。 常量和变量的定义按照论文中的谐音编写以便学者对照论文和模型进行学习。 因此对小白学习非常友好。 采用自己编写的龙格库塔算法对动力学方程求解。 输出了不同周期下内圈的质心运动轨迹内圈接触力等动态响应。轴承动力学模型复现这事儿说难不难但细节特磨人。咱先看核心方程内圈质心运动方程长得像俩弹簧拴着铁球左边是惯性项右边是接触力加阻尼。重点在于那个非线性赫兹接触力Qc的计算论文里写成QcKc*δ^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[i1] - 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[i1] 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转速转换容易翻车记住1rpm2π/60 rad/s。有个硕士生把转速直接当角频率用结果轨迹图变成心电图当场自闭。结果可视化时得注意采样间隔尤其观察周期性plt.figure(figsize(12,6)) plt.plot(t[::10], sol[::10,0], b, alpha0.7, labelX位移) # 降采样显示 plt.plot(t[::10], sol[::10,2], r, alpha0.7, labelY位移) plt.xlabel(时间(s)) plt.ylabel(位移(m)) plt.title(内圈质心运动轨迹) plt.legend()这里用了::10进行降采样不然20000个数据点直接绘图会卡成PPT。alpha参数调节透明度能更好显示轨迹重叠部分。最后说个实战经验当接触力出现负值时多半是赫兹接触力的符号处理有误。这时候应该检查delta计算是否带绝对值、sign函数应用是否正确。另外建议用scipy的solve_ivp做结果对比避免算法实现错误导致整组数据报废。
