简介基于ISO16281标准的圆柱滚子轴承刚度计算Python程序面向机械设计、轴承分析及转子动力学方向的工程师与研究者。程序将标准中的载荷-变形关系转化为可执行脚本输入基本几何与工况参数即可快速获得刚度结果适合用于轴承选型、转子系统建模或教学演示。压缩包内共1个文件为单个.py脚本压缩后仅1KB轻量精简、无额外依赖方便直接查看和修改算法逻辑。借助Python的数值处理能力该脚本能够替代传统手算流程减少公式推导与查表的繁琐同时保证计算过程透明、参数可调。目前已有222人学习/下载说明该程序在相关领域具备一定的参考价值通过该脚本读者可以结合ISO16281条文理解圆柱滚子轴承刚度的计算流程并以此为基础扩展不同工况下的参数化分析或嵌入其他仿真项目中作为计算模块使用。1. 为什么用 Python 解圆柱滚子轴承刚度从一个压缩包说起你拿到的文件叫 XYZhylinderBearing.rar这类压缩包打开后通常有两样东西一个装着滚子或滚道表面采样点的 .xyz 坐标文件一份计算轴承刚度的 Python 脚本。圆柱滚子轴承的刚度不是查表得到的固定值它取决于当前载荷、滚子数目、有效接触长度和轴承游隙同一套轴承在不同工况下径向刚度可能差 20% 以上。Python 适合做这件事是因为从散点数据到 Hertz 接触迭代再到刚度矩阵输出全流程都能用几十行代码闭环改一个参数就能重算一组曲线。这篇文章面向想把“轴承刚度”从经验值换成可计算量的工程师也面向刚在 VSCode、PyCharm 或 Linux 系统里搭好 Python 环境、准备跑通第一个算例的人。2. 圆柱滚子轴承刚度的计算模型Hertz 接触与滚子载荷分布2.1 线接触的弹性趋近量刚度为什么随载荷变圆柱滚子轴承里滚子与内圈、外圈滚道的接触属于线接触。Hertz 接触理论解决的是接触应力分布但工程算刚度时更常用的是载荷 Q 与弹性趋近量 δ 之间的幂函数关系。钢制圆柱滚子轴承常用的 Palmgren 形式写成δ C_delta * Q^n接触指数 n 在 0.9 附近C_delta 是跟滚子有效长度、材料弹性模量、滚子直径有关的系数。这里的单位必须自洽最常见的坑是把 Q 从 kgf 换成 N 后系数没有跟着换算结果位移差一个数量级。对 Q 求导立刻得到单对接触的切向刚度k dQ/dδ。因为有 0.9 这个小数指数k 并不恒定而是随 δ 增大慢慢变大。实际轴承载荷从空载增加到额定载荷时径向刚度通常会有明显爬升。因此任何把“轴承刚度”当成单一定值带入转子模型的做法都会在轻载工况低估、在重载工况高估支撑效果。# 单滚子接触刚度对载荷-位移的幂函数做数值微分 def roller_contact_stiffness(delta, C_delta, n_exp0.9): # delta: 滚子与滚道总弹性趋近量(mm)C_delta: 位移-载荷系数 h max(abs(delta) * 1e-6, 1e-10) q_plus ((delta h) / C_delta) ** (1.0 / n_exp) q_minus ((delta - h) / C_delta) ** (1.0 / n_exp) return (q_plus - q_minus) / (2.0 * h)代码里用中心差分避开了解析求导的指数运算h 跟随 delta 缩放避免小位移时出现除零。C_delta 和 n_exp 作为参数传入是为了让你换材质、换滚子修形时不用改函数本身。参数符号典型值说明载荷指数n_exp0.9Palmgren 线接触经验指数位移-载荷系数C_delta3.84e-5必须与力、位移单位制一致有效接触长度L_eff0.80.9×滚子全长端部修形后打折2.2 位移协调只有一部分滚子在承载单个滚子的刚度再高也要放回整圈滚子里面看。外圈固定、内圈受到径向力 Fr 时假设内圈沿载荷方向移动 δr那么第 j 个滚子的接触变形不是直接等于 δr而是其投影δ_j δr * cos(ψ_j) - Pd/2ψ_j 是这个滚子与载荷方向的夹角Pd 是轴承直径游隙。cos(ψ_j) ≤ 0 的滚子在受载侧对面理论上不承载δ_j ≤ 0 的滚子同样退出接触。只有载荷方向附近一组滚子真正受力这就是为什么圆柱滚子轴承的径向刚度比同尺寸球轴承高却对游隙更敏感。实际计算时不能给所有滚子平均分配载荷。给定一个试算的 δr每个承载滚子的 Q_j 由幂函数反解出来再校验力的平衡Fr - Σ Q_j * cos(ψ_j) 0这个单变量方程用 scipy.optimize.brentq 解是最省事的不用自己写 Newton-Raphson 的雅可比。下面这段就是最小可运行的内圈位移求解import numpy as np from scipy.optimize import brentq def solve_inner_ring_displacement(Fr, psi, Pd, Z, C_delta, n_exp0.9): # psi: 各滚子相对载荷方向的角度(rad)Pd: 直径游隙(mm) # 返回使力平衡的内圈径向位移(mm) def equilibrium(dr): delta_j dr * np.cos(psi) - Pd / 2.0 delta_j np.where(delta_j 0.0, delta_j, 0.0) q_j np.where(delta_j 0.0, (delta_j / C_delta) ** (1.0 / n_exp), 0.0) return float(np.sum(q_j * np.cos(psi)) - Fr) # 力的残差 return brentq(equilibrium, 1e-6, 2.0, xtol1e-12)brentq 要求两端点异号所以位移上界要先给足够大一般滚子轴承弹性变形在微米到几十微米量级2 mm 已经足够保守。xtol 设到 1e-12 是为了后面数值求导不引入迭代误差。2.3 从单滚子刚度到整机刚度矩阵求出一个载荷下的 δr 后把平衡条件对 δr 求导就是径向刚度 Krr。轴向刚度、角刚度的原理也一样只是位移方向换成轴向或转角再列相应方向的力或力矩平衡方程。若转子模型需要完整支撑信息常见做法是组装 5×5 对称矩阵把径向、轴向、角向耦合项一起输出只关心临界转速时取 Krr 和 Kaa 两个对角项就够。提示Palmgren 系数 C_delta 在不同手册里因单位制不同差很多。代码里先固定一套单位用简单拉伸试验或有限元算一个已知点反标 C_delta比直接抄公式常数更可靠。3. 用 Python 把 XYZ 数据变成圆柱滚子轴承刚度核心代码3.1 用 Pandas 读入滚子坐标并提取几何参数解压 XYZhylinderBearing.rar 后如果拿到的是 .xyz 散点文件第一步不是算刚度而是确认坐标系的轴向。圆柱滚子轴承建模时通常把旋转轴设为 z 轴x、y 是横截面坐标。先读进来再用柱坐标变换提取半径和角度这是最快的一次“体检”。import pandas as pd import numpy as np df pd.read_csv(roller.xyz, sepr\s, headerNone, names[x, y, z]) r np.hypot(df[x], df[y]) # 到轴线的半径 theta np.arctan2(df[y], df[x]) # 周向角度 R_surface float(np.median(r)) # 表面采样半径 L_eff float(df[z].max() - df[z].min()) # 有效接触长度(还需扣除倒角) L_eff * 0.85 # 考虑滚子端部修形的折算用 np.median 而不是 np.mean是为了避免个别离群点把半径估计带偏。L_eff 直接取 z 范围只能得到滚子全长真实接触长度要按滚子母线修形量打折常见做法是先乘 0.8~0.9 再通过刚度校核微调。这段代码不依赖特殊库pandas、numpy 在新装的 Python 环境里就能跑。我自己处理这类点云时第一步永远是按 z 分层、每层求半径再取中位数避免把倒角点和修形点混进有效直径。如果你现在还在为 Python 环境安装烦恼我的顺序是在 Linux 系统安装 Python 后立刻用 python3 -m venv .venv 建虚拟环境VSCode 配置 Python 环境时选这个解释器路径PyCharm 配置 Python 环境时也指向同一份。三件套 numpy/scipy/pandas 用 pip 装后面所有算例都是这套环境。文件类型常见列含义在刚度计算里的用途roller.xyzx y z 点坐标求滚子半径与有效长度inner.xyzx y z 点坐标内圈滚道半径算配合过盈race.xyzx y z 点坐标外圈滚道半径校核载荷分布3.2 用 brentq 求解内圈径向位移平衡几何参数有了接下来把上一章的平衡方程套上去。需要额外准备的是滚子位置角 psi如果压缩包里给出了装配图或滚子中心坐标用反正切算没有的话就假设滚子等距分布psi 2π * j / Z psi0psi0 由载荷方向与最近滚子的相对位置决定。Z 18 # 滚子数量 psi0 0.0 # 第一个滚子的初始角可调 psi psi0 2.0 * np.pi * np.arange(Z) / Z Fr 5000.0 # 工作径向载荷(N) Pd 0.012 # 直径游隙(mm)C3 游隙一般更大 C_delta 3.84e-5 # 按你采用的手册单位体系标定 dr_sol solve_inner_ring_displacement(Fr, psi, Pd, Z, C_delta) print(f内圈径向位移: {dr_sol * 1e3:.2f} um)psi0 这个参数很少被关注它决定了载荷方向是正对一个滚子还是正对两个滚子之间的空隙刚度结果因此会有小幅度周期变化对应实际轴承中的“滚动体通过频率”成分。如果做转子瞬态分析psi0 应该随转角更新而不是固定。这一段没有新函数但参数选择会直接影响结果。Pd 取 0.012 mm 代表普通组游隙换成负的 -0.005 mm 就变成预紧状态求解区域和承载滚子数都会改变。后面第 4 章专门讲游隙怎么改。3.3 中心差分与样条拟合把刚度曲线算稳单点位移不是最终目标刚度需要 dFr/dδr。最直接的办法是对一系列载荷分别求平衡位移再用 np.gradient 求导。这个办法的坑在于载荷取点太密差分放大迭代误差太稀曲线不够平滑。def bearing_stiffness_curve(loads, psi, Pd, Z, C_delta, n_exp0.9): drs [solve_inner_ring_displacement(F, psi, Pd, Z, C_delta, n_exp) for F in loads] drs np.array(drs) K_rr np.gradient(loads, drs) # dFr/dδr单位 N/mm return drs, K_rr loads np.geomspace(500, 50000, 12) drs, K_rr bearing_stiffness_curve(loads, psi, Pd, Z, C_delta) for d, K in zip(drs * 1e3, K_rr): print(f位移 {d:.2f} um径向刚度 {K:.1f} N/mm)np.gradient 默认对非均匀网格也做中心差分所以载荷用 geomspace 没问题。如果输出曲线有锯齿先别急着滤波检查平衡求解的 xtol再把载荷点改成 16~24 个多数锯齿来自载荷间隔跨度过大。到这一步你已经能从一包散点坐标走到一条 K-Fr 曲线。接下来能否用于工程取决于参数怎么设。4. 圆柱滚子轴承刚度实战参数游隙、预紧和转速怎么改4.1 直径游隙与负游隙的公式修正前面公式里的 Pd/2 是正游隙的处理方式即内圈先空走一段接触后才产生弹性变形。负游隙预紧时所有滚子都处于受压状态位移协调式要改写成δ_j δr * cos(ψ_j) δ_preδ_pre 是装配产生的初始压缩量。代码上的差别是把- Pd/2换成 preload_displacement同时去掉 np.where 里的零截断。这一改承载滚子数从“半圈”变成“全圈”刚度会明显抬升。很多轴承座振动问题就是预紧量被安装误差改掉后支撑刚度突变引起的。def solve_inner_ring_displacement_preload(Fr, psi, pre_delta, Z, C_delta, n_exp0.9): def equilibrium(dr): delta_j dr * np.cos(psi) pre_delta # 负游隙全部滚子接触 q_j (delta_j / C_delta) ** (1.0 / n_exp) return float(np.sum(q_j * np.cos(psi)) - Fr) return brentq(equilibrium, -1.0, 2.0) # 注意下界允许为负这里上下界区间要扩大到 -1.0因为预紧时内圈可以向载荷反方向移动一段。对预紧轴承谈“径向位移为 0”是不成立的它在无外载时已经有初始变形。4.2 高速工况下离心力对刚度的影响转速上去后滚子离心力会把滚子压向外圈内圈载荷变小、外圈载荷变大载荷分布不再关于载荷方向对称。离心力按 F_c m_c * ω² * D_m / 2 估算m_c 是单个滚子质量D_m 是节圆直径ω 是转动角速度。工程上常用等效径向力修正把每个滚子的离心力加到外圈接触方程里内圈仍按静态平衡算内外圈各自的刚度分开输出。这个效应在低转速、轻载时几乎可以忽略但在 d_m * n 超过 50 万d_m 是节圆直径 mmn 是转速 r/min时不能忽略。如果你的工况是高速主轴还在用静态刚度矩阵临界转速会算得偏保守或偏冒险取决于载荷方向。4.3 参数对照表与结果校验给出一组可直接起步的推荐参数参数符号初值说明滚子数Z18设计图或实测有效接触长度L_eff0.85×滚子全长端部修形后打折直径游隙Pd0.012 mm普通组C3 约 0.02~0.03预紧位移δ_pre0 或 0.005 mm仅负游隙时启用载荷指数n_exp0.9钢制滚子常规值载荷系数C_delta3.84e-5必须与单位制一致改参数后怎么判断结果真的对最有效的校验是渐进一致性把 Z 改大成 30、把 C_delta 改小 20%刚度曲线应平滑变化不能出现某个载荷点突变。更严一点和同一轴承的有限元静力分析对比误差在 5% 以内说明游隙和接触长度取得到位误差指向同一个方向则通常是 C_delta 或 L_eff 的系统偏差。载荷(N)本文算法位移(μm)有限元位移(μm)偏差1000待填待填5%5000待填待填5%20000待填待填5%注意brentq 求解时如果上下界都取正值预紧工况会直接报错。出现“f(a) 和 f(b) 必须异号”时先检查是否把负游隙当成正游隙处理了。5. 把算出的轴承刚度接进转子动力学模型5.1 用幂律拟合检查 K-Fr 曲线是否自洽理论上径向刚度应近似满足 K a * Fr^bb 约在 0.1。用 scipy.optimize.curve_fit 去拟合前面算出的曲线如果 b 偏离 0.1 太远多半是游隙或预紧的符号写反了而不是理论有问题。from scipy.optimize import curve_fit def power_law(F, a, b): return a * F ** b popt, _ curve_fit(power_law, loads, K_rr, p0[1e4, 0.1]) print(f拟合指数 b {popt[1]:.3f})这个技巧能在没有参考解的情况下快速筛掉低级错误。b 落在 0.05~0.15 之间时曲线形态可信b 落在 0.5 以上就去查载荷分布里的 psi 是否写成了角度制。5.2 导出 5×5 刚度矩阵给转子模型调用最后把 K 存成矩阵或轻量 CSV转子动力学软件里常见格式是 5 个自由度径向 x、径向 y、轴向 z、绕 x、绕 y。圆柱滚子轴承的交叉耦合项通常很小但预紧后不能直接设为 0。一次完整计算可以同时输出两个径向方向刚度再按工况写入文件供后续临界转速或不平衡响应分析使用。建议在 Linux 系统安装 Python 后把整套脚本放进 venv用 pip freeze 锁定依赖版本。VSCode 配置 Python 环境时加一个 Jupyter 扩展方便逐段看位移-载荷迭代过程。整个流程里最容易被忽略的不是 Hertz 公式而是坐标方向定义z 轴一旦反了psi 阵列全错结果却在数值上看起来合理。把上面两段组合成一个 check_curve.py白天算曲线、晚上批量跑参数扫描日志里记录每一组的 psi0 和游隙基本就能覆盖大多数轴承刚度复查需求。本文还有配套的精品资源点击获取