滚动轴承非线性动力学复现：从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_x0)、(W_y-mg)负号是重力方向(e) 是质量偏心距(\Omega) 是转轴角速度。我复现时用的是一组非常接近6205深沟球轴承的参数滚珠数 (N_b9)滚珠直径 (d7.94,\text{mm})节圆直径 (D39.04,\text{mm})径向游隙 (c_010,\mu\text{m})等效质量 (m0.6,\text{kg})接触刚度取 (K_b7.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 j1,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_{j1}^{N_b} \delta_j^{3/2} \cdot H(\delta_j) \cdot \cos\theta_j ][ F_y K_b \sum_{j1}^{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不是情怀而是性价比之选。方案局部截断误差每个步长计算导数成本适用场景显式Euler(O(h^2))1次演示教学实际工程极少用RK4(O(h^5))4次中高精度、非线性不极端生硬的系统solve_ivpRK45自适应自适应理论高次数不固定内部有误差控制黑盒快速计算但中间过程不可控轴承这个系统每个导数函数里要循环计算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}_{n1} \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})也就是微秒级步长。第二步做收敛性验证。分别用 (h2\times10^{-6})、(h1\times10^{-6})、(h5\times10^{-7}) 跑同一段稳态工况对比x-y轨迹的最大偏差。如果结果基本重合说明步长已经足够小。我实测中(h1\times10^{-6}) 和 (h5\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[i1] rk4_step(derivatives, t_i, s_arr[i], h) t_arr[i1] (i 1) * h这个主循环里有一个很多人都会犯的隐性错误把积分初始阶段的瞬态响应直接拿去画轨迹。对于强非线性系统如果初值设置不当前0.05~0.1秒内系统会有一段剧烈的过渡过程轨迹图上会多出一大团乱线。我的做法是总时长至少积分0.3秒以上画图时只取最后0.1秒的数据。初值方面我直接用零点起算(x0, y0, v_x0, v_y0)。但要注意零点起算后系统会在重力下迅速下沉进入接触状态这个过程正是瞬态响应的一部分。如果想更快进入稳态可以把初值设为静态平衡点附近但论文复现时我建议还是从零点起算减少人为干预。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, lw0.6, colorb) plt.xlabel(xc (m)) plt.ylabel(yc (m)) plt.axis(equal) plt.title(Mass Center Orbit) plt.grid(True, alpha0.3)plt.axis(equal)这行一定要加。不加的话matplotlib会自动拉伸坐标轴原本接近圆形的轨迹会被拉成椭圆肉眼判断形状时容易产生严重误导。4.2 等比例坐标与轨迹后处理画轨迹时还有几个实操细节值得说。第一轨迹线宽要调低。因为稳态轨迹是几十万个点首尾相连线宽太大时会糊成一团黑色看不到任何细节。我通常用lw0.4~0.8。第二可以用颜色映射转角度。按时间着色用scatter函数把每个点按时间顺序从蓝色渐变到红色能直观看出转子是正进动还是反进动以及运动是否存在漂移。这个信息在轴承故障诊断里很有价值。第三看轨迹的均值中心。由于重力让转子在垂直方向有一个静态偏移轨迹的几何中心通常不是原点。标准圆轨迹说明系统接近各向同性如果轨迹明显扁平说明两个方向的等效刚度差异很大很可能是轴承预紧或游隙设置带来的。4.3 频域验证轨迹对不对看频谱就知道画完轨迹我还会做一个频域验证这一步能快速判断复现是否成功。对稳态的 (x(t)) 序列做FFTxf np.fft.rfft(x) freq np.fft.rfftfreq(len(x), dh) 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_r50,\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}}) 到底取什么有人取系统一阶固有频率有人直接取转频。不同取法对轨迹形态影响很大。复现时我建议优先参考原文给出的阻尼设置如果原文没给就先用一个中等阻尼比 (\zeta0.01\sim0.05)然后做敏感性分析——把阻尼比从0.01扫到0.1观察轨迹从“瘦”到“胖”的变化趋势。这个趋势本身也是有价值的分析内容能写进复现报告里。5.4 游隙和过盈(\delta_j) 的符号约定游隙 (c_0) 的定义域正负在不同论文里可能完全相反。我们这里用的约定是(\delta_j 0) 表示滚动体受压变形。若轴承处于预紧过盈状态游隙为负值模型里的 (c_0) 要写成负数或者用另一个变量表示过盈量。这个符号一旦搞反(F_x)、(F_y) 的方向不会错但接触力的出现时机全错时变刚度特性就完全变味了。一个快速判断方法给系统加一个很小的静态力观察稳态位移方向。如果位移方向和力的方向相反基本可以确定游隙符号写反了。5.5 判断复现成功的量化指标复现不是“画出来像就行”。真正说服自己的指标包括以下几项指标合理范围/特征稳态轨迹形状接近准椭圆或多瓣花瓣形中心偏移约 (mg/K_{\text{avg}}) 量级位移幅值量级十几到几十微米视转速和不平衡量而定FFT主频转频处有峰BPFO附近有边频调制步长收敛性(h) 减半后轨迹最大偏差小于0.1(\mu\text{m})质心轨迹与轴心轨迹差异差异约为偏心距 (e) 对应的旋转圆复现成功后我一般会把稳态轨迹、时间历程、FFT谱图三张图放在一起作为完整交付。这样无论在论文附录还是技术报告里都有充分证据支持你的结论。说回我开头遇到的那个NaN问题——最后查明就是 (K_b) 单位搞错导致接触力远小于正常值系统在游隙里自由漂移位移发散。修正后立刻稳定。复现论文这件事技术门槛其实不高真正难的是对每个参数、每个符号都保持“怀疑并验证”的态度。希望这篇踩坑记录能帮你少走几步弯路。

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

相关推荐

ToF深度相机实战：从硬件到应用的全链路解析

Pyright 导入解析与打包机制深度解析：Import Resolution and Packaging 全指南

Data-Science-For-Beginners 分布可视化课后实战：用直方图与密度图为新数据集编写数据故事 Notebook

更多资讯

PeakTech P1260台式示波器：12位ADC+触摸屏的产线级实用主义选择

SpringBoot实现企业级Wiki系统的RBAC权限管理

FOC驱动小体积高扭矩瓶颈：MCU时间精度与MOS开关损耗硬约束

VS Code搭建STM32开发环境：从安装到编译烧录全流程

汽车电子精密制造数字化转型：ONES解决方案解析

Python电商推荐系统实战：从算法到毕业设计

今日更新

51单片机开环控制磁阻传感器的硬件匹配与代码实现

AI Agent安全操作数据库：SQL注入防御与参数化查询实战

Android逆向实战：还原激励类App签到接口的签名算法

本周精选

AI驱动的自动化专利分析平台技术解析

LeetCode-Go 题解：1252. Cells with Odd Values in a Matrix（奇数值单元格计数）双解法剖析

在线答疑系统Java毕设实战：Spring Boot前后端分离与状态流转

本月热门

基于YOLO的焊缝缺陷检测：从数据标注到工业部署全流程实践

分布式服务的决策记录

创业产品的首版范围

读完文章，想聊聊您的网站？