滚动轴承非线性动力学复现：从RK4算法到质心运动轨迹的完整实战

复现一篇轴承动力学方向的论文，最难受的不是公式看不懂，而是公式写得明明白白，跑出来结果死活对不上。前阵子我复现了一篇滚动轴承非线性动力学分析的paper，目标很明确：输出转子质心运动轨迹，验证论文里的轴心轨迹图。结果数值一度直接飞成NaN，排查到最后才发现是接触刚度 (K_b) 的单位换算差了4个数量级。这种坑，论文里不会写，审稿人也不会问，只有自己动手敲代码时才会撞上。

这篇文章就把整个复现过程摊开来讲：从轴承动力学模型的方程怎么落到可编码的形式，到四阶龙格库塔算法（RK4）的选型理由和步进公式，再到质心运动轨迹的绘图和频域验证，最后集中盘点几个真正折磨人的细节。适合正在做转子动力学、滚动轴承故障诊断方向的研究生，也适合想理解轴承非线性振动仿真逻辑的工程师。我不打算只给代码，更想讲清楚每步选择背后的原因——因为复现论文这件事，参数一致才是最终的评判标准。

1. 复现前先把方程按“可编码”的标准重写一遍

很多论文的公式排版很漂亮，下标和角标满天飞，真到写代码时你会发现，有一半的符号是模板生成的冗余，另一半则有歧义。所以我的第一步不是急着写RK4，而是把模型方程重新整理成一份编码能直接对照的清单。

1.1 二自由度集中质量模型的物理图像

这类轴承动力学论文里最常见的，是一个二自由度的集中质量模型。它把滚动轴承-转子系统简化成一个等效质量 (m)，这个质量在轴承平面内沿水平方向 (x) 和垂直方向 (y) 运动。方程长这样：

[ m\ddot{x} + c\dot{x} + F_x = W_x + m e \Omega^2 \cos(\Omega t) ]

[ m\ddot{y} + c\dot{y} + F_y = W_y + m e \Omega^2 \sin(\Omega t) ]

其中 (c) 是系统阻尼，(F_x) 和 (F_y) 是滚动体与滚道接触产生的Hertz接触力在水平和垂直方向的合力，(W_x) 和 (W_y) 是静载荷分量（轴水平放置时通常 (W_x=0)、(W_y=-mg)，负号是重力方向），(e) 是质量偏心距，(\Omega) 是转轴角速度。

我复现时用的是一组非常接近6205深沟球轴承的参数：滚珠数 (N_b=9)，滚珠直径 (d=7.94,\text{mm})，节圆直径 (D=39.04,\text{mm})，径向游隙 (c_0=10,\mu\text{m})，等效质量 (m=0.6,\text{kg})，接触刚度取 (K_b=7.055\times10^9,\text{N/m}^{3/2})。这些参数在相关文献中很常见，但要注意——每篇论文对“等效质量”的定义不一样，有的含转子分配质量，有的只算轴承座局部质量，复现时需要根据论文场景判断。

1.2 Hertz接触力与滚动体间歇接触

滚动轴承最关键的非线性源，来自滚动体与滚道的Hertz接触。每个滚动体的角位置随时间变化：

[ \theta_j = \frac{2\pi(j-1)}{N_b} + \omega_c t,\quad j=1,2,\dots,N_b ]

这里 (\omega_c) 是保持架角速度，外圈固定时近似为：

[ \omega_c \approx \frac{\Omega}{2}\left(1-\frac{d}{D}\right) ]

第 (j) 个滚动体在径向方向上的接触变形为：

[ \delta_j = x\cos\theta_j + y\sin\theta_j - c_0 ]

注意这里有个极其容易忽略的约束条件：Hertz接触力只在变形为正时存在。如果 (\delta_j \le 0)，滚动体与滚道脱离接触，该滚动体不贡献任何力。所以实际编码时要加一个Heaviside判断。接触力分量写为：

[ F_x = K_b \sum_{j=1}^{N_b} \delta_j^{3/2} \cdot H(\delta_j) \cdot \cos\theta_j ]

[ F_y = K_b \sum_{j=1}^{N_b} \delta_j^{3/2} \cdot H(\delta_j) \cdot \sin\theta_j ]

这个表达式本身不复杂，但它带来的系统行为非常“野”：每个滚动体进入承载区、离开承载区的瞬间，等效刚度都会突变。转起来之后，系统受到的激励不是单一频率的正弦力，而是包含转频、滚动体通过频率及其倍频的周期时变刚度。这也直接决定了后面的数值积分为什么不能随便选个Euler法就上。

1.3 把二阶方程组改写成状态空间形式

RK4本质上只能处理一阶微分方程组，所以我们必须先把两个二阶方程改写成一阶状态空间。取状态向量：

[ \mathbf{s} = [x,\ y,\ \dot{x},\ \dot{y}]^T ]

那么一阶导数为：

[ \dot{\mathbf{s}} = \begin{bmatrix} \dot{x} \ \dot{y} \ \ddot{x} \ \ddot{y} \end{bmatrix}

\begin{bmatrix} v_x \ v_y \ \dfrac{1}{m}(W_x + m e \Omega^2 \cos\Omega t - c v_x - F_x) \ \dfrac{1}{m}(W_y + m e \Omega^2 \sin\Omega t - c v_y - F_y) \end{bmatrix} ]

这一步看起来只是数学上的搬运，但它是编码前最重要的“翻译”动作。很多新手直接对每个加速度分量分别写函数，结果状态向量顺序乱掉，RK4里 (k_2)、(k_3) 更新时就全错了。我的建议是：先把状态向量里每个分量的意义写在注释里，再写导数函数，保持顺序完全一致。

2. 四阶龙格库塔为什么是这篇paper复现的正确答案

标题里带着“龙格库塔算法”，说明这套数值积分方案本身就是复现核心。但为什么偏偏是RK4？我一开始也想省事，直接用Python的scipy.integrate.solve_ivp，结果并不理想——不是说不能算，而是论文复现时你很难解释清楚“容差设置、步长自适应策略”对结果到底有多大影响。自己手写RK4，所有参数都在掌控之中，复现文献结果时更有说服力。

2.1 三个候选方案：Euler、ode45、RK4

先做个直观对比，让大家明白RK4不是情怀，而是性价比之选。

轴承这个系统，每个导数函数里要循环计算9个滚动体的接触力，本身就不便宜。Euler要保证精度必须把步长压到极小，计算量反而更大。RK4在每个步长内把系统状态推进得足够准，四步加权平均的误差特性非常适合这种存在高频接触切换的振荡系统。

更重要的是，论文复现时需要把“数值方法”当作报告的一部分写清楚。如果审稿人或导师问你“用了什么积分器、步长多少、截断误差多少”，手写RK4能直接答得明明白白，丢一个solve_ivp回去反而显得像黑盒实验。

2.2 RK4步进公式的记忆方式与局部截断误差

四阶龙格库塔的步进公式是：

[ k_1 = f(t_n, \mathbf{s}_n) ]

[ k_2 = f\left(t_n + \frac{h}{2},\ \mathbf{s}_n + \frac{h}{2}k_1\right) ]

[ k_3 = f\left(t_n + \frac{h}{2},\ \mathbf{s}_n + \frac{h}{2}k_2\right) ]

[ k_4 = f\left(t_n + h,\ \mathbf{s}_n + h k_3\right) ]

[ \mathbf{s}_{n+1} = \mathbf{s}_n + \frac{h}{6}(k_1 + 2k_2 + 2k_3 + k_4) ]

我记忆这个公式的方法很简单：一辆车的速度不是恒定的，想在时间 (h) 内准确估算位移，不能只看起点速度。最好在起点、两个中间点、终点各测一次速度，然后按“1、2、2、1”的权重求平均。(k_1) 是起点斜率，(k_2) 是用起点斜率试探到中点的修正斜率，(k_3) 是用修正斜率再测中点的更优斜率，(k_4) 是用更优斜率试探到终点的斜率。这么一来，RK4的每一步包含了四次导数求值，能捕捉到系统在中途的变化，而不是像Euler那样“一条道走到黑”。

这个权重公式对应的是Simpson数值积分思想，不需要死记硬背，自己可以现场推一遍。它的全局误差是 (O(h^4))，局部误差是 (O(h^5))。这也是“四阶”这个说法的来源。

值得注意的是，(k_2) 和 (k_3) 都算中点斜率，但一个基于 (k_1) 试探，一个基于 (k_2) 修正，两者物理意义不同，只有多次编码写错过“把 (k_3) 写成用 (k_1) 算”的人，才会明白为什么这一步顺序不能乱。

2.3 步长怎么拍：从特征频率到收敛性验证

步长选择是复现中最容易蒙圈的地方。我用的实用方法分两步。

第一步，估算系统里最高的相关频率。以6205轴承、3000rpm为例：

• 转频 (f_r = 50,\text{Hz})
• 保持架特征频率 (f_c \approx 19.9,\text{Hz})
• 外圈故障特征频率 (f_{BPFO} = N_b \times f_c \approx 179,\text{Hz})

但由于滚动体进出承载区会产生冲击，实际响应里往往有更高次的谐波分量，尤其接触刚度非线性强时，能量可以延伸到kHz级别。为了保证积分稳定，步长建议取最高关注频率周期的1/100到1/1000。针对这个系统，我选了 (h = 1\times10^{-6},\text{s})，也就是微秒级步长。

第二步，做收敛性验证。分别用 (h=2\times10^{-6})、(h=1\times10^{-6})、(h=5\times10^{-7}) 跑同一段稳态工况，对比x-y轨迹的最大偏差。如果结果基本重合，说明步长已经足够小。我实测中，(h=1\times10^{-6}) 和 (h=5\times10^{-7}) 的轨迹差异在微米量级以下，可以接受。这个收敛性验证一定要做，否则论文复现出来的轨迹可能只是数值误差的产物。

3. RK4编码实践：从derivatives到积分主循环

代码不长，但每一块都有讲究。我直接给出一个精简版，跑通之后你可以按需扩展。

3.1 参数区：全部使用SI单位

import numpy as np import matplotlib.pyplot as plt # ---------- 轴承与系统参数（SI单位） ---------- Nb = 9 # 滚动体数量 d = 7.94e-3 # 滚动体直径 m D = 39.04e-3 # 节圆直径 m c0 = 10e-6 # 径向游隙 m Kb = 7.055e9 # Hertz接触刚度 N/m^1.5 m = 0.6 # 等效集中质量 kg # 转速与激励 rpm = 3000 Omega = 2 * np.pi * rpm / 60 # 转轴角速度 rad/s fr = rpm / 60 # 转频 Hz omega_c = Omega / 2 * (1 - d / D) # 保持架角速度 rad/s ecc = 20e-6 # 质量偏心距 m g = 9.81 Wx = 0.0 Wy = -m * g # 重力垂直向下 # 阻尼，按线性阻尼处理 zeta = 0.02 c_damp = 2 * zeta * m * Omega # 简化估计，实际可按系统模态修正

单位统一是复现的第一道坎。把 (d) 写成 (7.94e-3)，而不是 (7.94)，是因为整个公式都以米为基准。后面 (K_b) 的单位是 (\text{N/m}^{3/2})，如果你的参考论文给的是 (\text{N/mm}^{3/2})，记得换算，具体怎么换我放到最后一节详细说。

3.2 导数函数里最容易算错的一项

def bearing_force(x, y, t): fx = 0.0 fy = 0.0 for j in range(Nb): theta = 2 * np.pi * j / Nb + omega_c * t delta = x * np.cos(theta) + y * np.sin(theta) - c0 if delta > 0.0: force = Kb * delta**1.5 fx += force * np.cos(theta) fy += force * np.sin(theta) return fx, fy def derivatives(t, s): x, y, vx, vy = s fx, fy = bearing_force(x, y, t) ax = (Wx + m * ecc * Omega**2 * np.cos(Omega * t) - c_damp * vx - fx) / m ay = (Wy + m * ecc * Omega**2 * np.sin(Omega * t) - c_damp * vy - fy) / m return np.array([vx, vy, ax, ay])

这段代码最容易错的地方不在RK4本身，而在Hertz接触力的分量映射。很多新手会把 (\delta_j^{3/2}) 当成标量力，然后直接加到全局力上。实际上每个滚动体的接触力方向沿径向，必须投影到 (x)、(y) 两个方向，才能和方程里的位移坐标对应。如果不投影，轨迹形状铁定不对。

另外一个细节是 (\delta_j) 的符号判断。if delta > 0.0这个判断就是前文说的Heaviside条件。实测中这个判断对性能影响很大，因为每个积分步长都要计算9次。如果未来想提速，可以用NumPy向量化一次算9个滚动体的合力，但初版复现最好保持循环结构，逻辑清楚比性能重要得多。

3.3 主循环与数据采集：别把瞬态段混进来

def rk4_step(f, t, s, h): k1 = f(t, s) k2 = f(t + 0.5*h, s + 0.5*h*k1) k3 = f(t + 0.5*h, s + 0.5*h*k2) k4 = f(t + h, s + h*k3) return s + (h/6.0)*(k1 + 2.0*k2 + 2.0*k3 + k4) # ---------- 积分主循环 ---------- t_total = 0.5 # 总时长 s h = 1e-6 # 步长 s n_steps = int(t_total / h) t_arr = np.zeros(n_steps + 1) s_arr = np.zeros((n_steps + 1, 4)) s_arr[0] = [0.0, 0.0, 0.0, 0.0] for i in range(n_steps): t_i = i * h s_arr[i+1] = rk4_step(derivatives, t_i, s_arr[i], h) t_arr[i+1] = (i + 1) * h

这个主循环里有一个很多人都会犯的隐性错误：把积分初始阶段的瞬态响应直接拿去画轨迹。对于强非线性系统，如果初值设置不当，前0.05~0.1秒内系统会有一段剧烈的过渡过程，轨迹图上会多出一大团乱线。我的做法是总时长至少积分0.3秒以上，画图时只取最后0.1秒的数据。

初值方面，我直接用零点起算：(x=0, y=0, v_x=0, v_y=0)。但要注意，零点起算后系统会在重力下迅速下沉，进入接触状态，这个过程正是瞬态响应的一部分。如果想更快进入稳态，可以把初值设为静态平衡点附近，但论文复现时我建议还是从零点起算，减少人为干预。

4. 质心运动轨迹输出：从状态向量到x-y相平面

积分完成后，状态向量每一行都记录了系统在不同时刻的 (x, y, v_x, v_y)。所谓质心运动轨迹，直观来说就是转子质心在轴承平面内的运动路线，画出来通常是一团类似椭圆或花瓣形的曲线。

4.1 什么是“质心运动轨迹”，它和轴心轨迹差了一个偏心矢量

这里必须澄清一个概念：多数文献说“轴心轨迹”时，指的是轴颈几何中心的位置变化，也就是我们积分得到的 (x(t)) 和 (y(t))。而标题里写的是“质心运动轨迹”，严格讲是转子质心 (O_c) 的运动轨迹。质心和几何中心之间差了一个偏心距矢量：

[ x_c(t) = x(t) + e\cos(\Omega t) ]

[ y_c(t) = y(t) + e\sin(\Omega t) ]

所以在画图时，我会同时输出两组轨迹：一组是几何中心轨迹 (x(t)-y(t))，另一组是基于偏心修正的质心轨迹 (x_c(t)-y_c(t))。两者形态相似，但质心轨迹多了一个以转频旋转的偏心圆分量，通常轨迹整体会显得更“胖”一些。复现论文时一定要先看清楚作者画的是哪一种，别拿轴心轨迹去硬比质心轨迹。

# ---------- 取稳态段 ---------- start_idx = int(0.3 / h) # 丢弃前0.3s瞬态 x = s_arr[start_idx:, 0] y = s_arr[start_idx:, 1] t_steady = t_arr[start_idx:] # 质心轨迹 xc = x + ecc * np.cos(Omega * t_steady) yc = y + ecc * np.sin(Omega * t_steady) # 绘制 plt.figure(figsize=(6, 6)) plt.plot(xc, yc, lw=0.6, color='b') plt.xlabel('xc (m)') plt.ylabel('yc (m)') plt.axis('equal') plt.title('Mass Center Orbit') plt.grid(True, alpha=0.3)

plt.axis('equal')这行一定要加。不加的话，matplotlib会自动拉伸坐标轴，原本接近圆形的轨迹会被拉成椭圆，肉眼判断形状时容易产生严重误导。

4.2 等比例坐标与轨迹后处理

画轨迹时还有几个实操细节值得说。

第一，轨迹线宽要调低。因为稳态轨迹是几十万个点首尾相连，线宽太大时会糊成一团黑色，看不到任何细节。我通常用lw=0.4~0.8。

第二，可以用颜色映射转角度。按时间着色，用scatter函数把每个点按时间顺序从蓝色渐变到红色，能直观看出转子是正进动还是反进动，以及运动是否存在漂移。这个信息在轴承故障诊断里很有价值。

第三，看轨迹的均值中心。由于重力让转子在垂直方向有一个静态偏移，轨迹的几何中心通常不是原点。标准圆轨迹说明系统接近各向同性；如果轨迹明显扁平，说明两个方向的等效刚度差异很大，很可能是轴承预紧或游隙设置带来的。

4.3 频域验证：轨迹对不对，看频谱就知道

画完轨迹，我还会做一个频域验证，这一步能快速判断复现是否成功。对稳态的 (x(t)) 序列做FFT：

xf = np.fft.rfft(x) freq = np.fft.rfftfreq(len(x), d=h) amp = np.abs(xf) / len(x) * 2 plt.figure(figsize=(8, 3)) plt.plot(freq, amp) plt.xlim(0, 1000) plt.xlabel('Frequency (Hz)') plt.ylabel('Amplitude')

对于健康的滚动轴承-转子系统，频谱中通常有两个显著特征：

• 转频 (f_r=50,\text{Hz}) 处有明显峰值，这是质量偏心激励直接产生的；
• 在 (f_{BPFO}\approx179,\text{Hz}) 附近出现调制边频带，这是滚动体通过时变接触刚度的结果。

如果复现结果里完全没有 (f_{BPFO}) 附近的分量，那基本可以断定模型或者参数有误。反过来，如果频谱被一大堆高频噪声淹没，那多半是步长太大、积分不稳定，需要减小步长后再试。

5. 复现中真正折磨人的五个细节

这一节的内容不是从哪篇paper里看来的，全是实际敲代码踩出来的经验。每一条都可能让结果从“漂亮轨迹”变成“一坨废线”。

5.1 接触刚度 (K_b) 单位换算：一错就是4个数量级

这是最坑的、没有之一。很多论文接触刚度给的是 (7.055\times10^7,\text{N/mm}^{3/2})，但我们的方程要用SI单位，需要换算。推导如下：

[ 1,\text{N/mm}^{3/2} = 1,\text{N}/(10^{-3},\text{m})^{3/2} = 10^{4.5},\text{N/m}^{3/2} \approx 31623,\text{N/m}^{3/2} ]

所以如果文献数值是 (7.055\times10^7,\text{N/mm}^{3/2})，换算成SI就是大约 (2.23\times10^{12},\text{N/m}^{3/2})。如果你忘了换算，直接拿 (7.055\times10^7) 当SI数值用，接触力会凭空缩小4个数量级，系统的静变形几乎为零，轨迹全部飘在游隙里，结果当然是一团乱。我在复现初期就栽过这个跟头，建议所有论文参数统统一到SI后再进代码，不要中途混合使用。

5.2 初值、舍弃周期数与前后处理

初值问题前面提过，再补充一个量化标准：判断系统是否进入稳态，可以观察最后两个转频周期的轨迹是否重合。方法是在绘图前，把稳态段按转周期切片，分别画在同一张图上，如果相邻两个周期的轨迹几乎重合，说明可以拿这段数据出图了。

同时注意积分时长的选择。积分太短，瞬态占比高；积分太长又浪费算力。我一般先试0.3秒，看轨迹稳定后截取后0.1秒，如果还不稳定就加长到0.5秒。另一个实用技巧是画x(t)的时间历程曲线，观察振幅是否在某个稳定水平附近波动，如果是，就说明系统已经进入稳态。

5.3 阻尼怎么给：线性阻尼近似的适用范围

轴承系统中的阻尼机制其实非常复杂，包括材料阻尼、油膜阻尼、接触迟滞等。绝大多数论文做非线性建模时都不细究阻尼机理，而是用一个等效线性阻尼 (c) 糊上去。我们的代码里用阻尼比估计：

[ c = 2\zeta m \Omega_{\text{ref}} ]

问题在于 (\Omega_{\text{ref}}) 到底取什么？有人取系统一阶固有频率，有人直接取转频。不同取法对轨迹形态影响很大。复现时我建议优先参考原文给出的阻尼设置；如果原文没给，就先用一个中等阻尼比 (\zeta=0.01\sim0.05)，然后做敏感性分析——把阻尼比从0.01扫到0.1，观察轨迹从“瘦”到“胖”的变化趋势。这个趋势本身也是有价值的分析内容，能写进复现报告里。

5.4 游隙和过盈：(\delta_j) 的符号约定

游隙 (c_0) 的定义域正负在不同论文里可能完全相反。我们这里用的约定是：(\delta_j > 0) 表示滚动体受压变形。若轴承处于预紧过盈状态，游隙为负值，模型里的 (c_0) 要写成负数，或者用另一个变量表示过盈量。这个符号一旦搞反，(F_x)、(F_y) 的方向不会错，但接触力的出现时机全错，时变刚度特性就完全变味了。

一个快速判断方法：给系统加一个很小的静态力，观察稳态位移方向。如果位移方向和力的方向相反，基本可以确定游隙符号写反了。

5.5 判断复现成功的量化指标

复现不是“画出来像就行”。真正说服自己的指标包括以下几项：

复现成功后，我一般会把稳态轨迹、时间历程、FFT谱图三张图放在一起，作为完整交付。这样无论在论文附录还是技术报告里，都有充分证据支持你的结论。

说回我开头遇到的那个NaN问题——最后查明就是 (K_b) 单位搞错导致接触力远小于正常值，系统在游隙里自由漂移，位移发散。修正后立刻稳定。复现论文这件事，技术门槛其实不高，真正难的是对每个参数、每个符号都保持“怀疑并验证”的态度。希望这篇踩坑记录能帮你少走几步弯路。