基于Harris理论的滚珠轴承准动力学MATLAB模型与载荷分布计算

简介：一套基于Harris理论的滚珠轴承准动力学模型MATLAB代码，面向机械、计算机、电子信息工程及数学专业的学生与研究者，适用于课程设计、期末大作业和毕业设计中的动力学仿真分析。代码兼容MATLAB 2014至2021a多个版本，采用参数化编程，参数可随时调整，注释详细，思路清晰，便于二次开发与学习。压缩包共13个文件，以9个.m脚本为主干，涵盖接触面积、油膜系数、椭圆求解等关键计算模块，另含3个备份文件和1个说明文档，包体仅12KB，轻量紧凑。已有37人浏览学习，附赠可直接运行的案例数据，无需额外收集即可复现结果。通过修改轴承参数，可直观观察不同工况下的动力学响应，既能加深对Harris理论准动力学模型的理解，也能提升MATLAB建模与调试能力，非常适合作为课程项目或毕业设计的参考资料。

1. 基于Harris理论的MATLAB准动力学模型：轴承刚度不再是黑匣子

用MATLAB算滚珠轴承的等效刚度，最怕的不是公式长，而是把Harris理论搬进准动力学模型时，几个量纲和迭代参数连续翻车。这份基于Harris理论的滚珠轴承准动力学模型MATLAB代码，把Hertz接触变形、钢球载荷分布和Newton-Raphson迭代串成一条完整链路，解决的是转子动力学分析里最现实的问题：轴承刚度矩阵和钢球载荷分布到底怎么算、怎么验。它不像有限元那样动辄几百万自由度，也不像纯经验公式那样查表查到怀疑人生，而是把Harris的滚动轴承力学分析经典理论落成了一套能跑、能改、能出图的脚本。适合做转子临界转速计算、轴承故障仿真、轴承选型校核的工程师，也适合正在啃滚动轴承拟静力学的研究生——拿到手不是背公式，而是直接改参数看结果。

2. Harris接触理论与准动力学平衡方程：核心公式与自由度取舍

2.1 Hertz点接触与Q-δ非线性关系

Harris理论的基础是Hertz弹性接触理论。滚珠与内外滚道之间是典型的点接触，接触区域的应力分布是椭球形的，而接触变形δ与载荷Q之间呈非线性关系：

Q = K × δ ^ 3/2

这个1.5次方是整个准动力学模型里最核心的非线性来源。K是载荷-变形常数，由接触体的材料弹性模量、泊松比和接触点的曲率半径共同决定。实际计算时K不是一个常数，它和接触点的主曲率和函数Σρ有关，曲率差F(ρ)决定Hertz接触椭圆的长短轴比，再由椭圆积分查表或数值求解得到K。

新手最容易出错的一点：K的单位是N/mm^1.5，不是N/mm。这意味着计算K时，弹性模量如果用GPa，必须先换算成MPa或N/mm²——1 GPa = 1000 N/mm²，漏掉这一步K会整整差三个数量级。我见过不少跑出来的载荷分布图完全违背常识，最后查到根因就是E的单位没换。

工程上还有一种常见做法是查Harris书里的K系数表，但表格只覆盖特定材料和曲率组合，换轴承型号就不准了。代码里直接按材料参数和几何参数现场算K，通用性更好，换轴承尺寸只需要改曲率半径数组。

2.2 从静力学到准动力学：五自由度模型的建立

准动力学模型的全称是quasi-static model，介于纯静力学和完全动力学之间。它忽略保持架碰撞、润滑剂瞬态阻尼这些高频效应，但把每个钢球的接触变形和载荷分布通过平衡方程严格求解。对大多数转子动力学分析来说，这种折中精度已经足够。

坐标系取轴承中心为原点，五个自由度分别是：两个径向平动x、y，一个轴向平动z，两个角向摆动θx、θy。五个自由度对应的外部载荷是Fx、Fy、Fz、Mx、My。每个钢球的位置由方位角ψj确定：

ψj = 2π × (j-1) / Z + ψ0

这里Z是钢球数，ψ0是第一个钢球的初始方位角，通常取0（即3点钟方向）。外圈固定时，钢球随保持架以约0.5倍轴转速公转，但在准动力学模型里，只要外载荷不随时间快速变化，公转带来的离心力和陀螺力矩可以折算成附加力，不需要做瞬态时域积分。

每个钢球的接触变形δj通过几何映射得到。以深沟球轴承为例，内圈在x、y、z方向平移后，第j个钢球处的法向趋近量是：

δj = x × cos(ψj) + y × sin(ψj) + z - （径向游隙影响项）+ 角向摆动修正项

注意这里径向游隙不能简单减一个常数。Harris的处理方法是：只有当δj > 0时钢球才受载，δj ≤ 0时该钢球载荷为0。这直接导致载荷分布在一个有限的角度区间内，而不是整个圆周都有钢球受力。载荷分布角ψ_l就是由游隙和径向位移共同决定的，这是一个隐式关系，必须放进迭代里一起解。

2.3 接触角变化：角接触球轴承的隐式耦合

深沟球轴承的接触角基本不变，但角接触球轴承不同：轴向预紧会使接触角从出厂名义值α0增大到工作接触角α。而接触角本身又取决于轴向位移和径向位移的组合，于是出现了一个隐式耦合——位移决定接触角，接触角又反过来决定位移，必须迭代求解。

代码里的处理思路是把每个钢球的接触角αj也作为局部变量，在每个迭代步内单独更新：先根据当前位移算几何关系，更新αj，再算δj，然后合成Qj，最后叠加成五个方向的合力与合力矩，和外部载荷做差得到残差。这个嵌套结构是准动力学模型代码最微妙的地方——主迭代是位移的Newton-Raphson，内层是每个钢球接触角的局部更新，两层的收敛阈值要区分开，否则会出现外层已经收敛但内层还在振荡的情况。

3. MATLAB代码结构与主流程：从输入参数到载荷分布

3.1 文件组织与运行入口

拿到这份代码后，先看目录结构再动手跑。完整的准动力学求解器通常按这样的模块划分，这是滚动轴承数值计算的常见工程做法：

bearing_harrris_qs/ ├── main.m # 主脚本，定义工况并调用求解器 ├── params_bearing.m # 轴承几何与材料参数 ├── compute_K.m # 计算Hertz接触刚度常数K ├── compute_delta.m # 由位移向量求每个钢球的接触变形 ├── compute_residual.m # 计算五个方向的力/力矩残差 ├── nr_solve.m # Newton-Raphson主迭代 └── plot_results.m # 后处理：载荷分布图、刚度曲线

main.m是唯一需要手动运行的脚本，运行前按轴承型号和工况修改params_bearing.m里的参数即可。输出包括每个钢球的法向载荷Qj、工作接触角、五个自由度的位移解，以及重构出的刚度矩阵。

3.2 输入参数定义与单位约定

params_bearing.m里的参数块长这样：

% ========== 轴承几何参数 ========== Z = 12; % 钢球数 db = 8.731; % 钢球直径，mm dm = 38.5; % 节圆直径，mm ri = 4.762; % 内圈沟曲率半径，mm re = 4.762; % 外圈沟曲率半径，mm alpha0 = 15.0; % 名义接触角，deg，角接触球轴承取15/25/40 % ========== 材料参数（注意单位） ========== E = 207000; % 滚动体和滚道弹性模量，N/mm^2（即MPa） nu = 0.3; % 泊松比 % 两个接触体材料相同，若不同分别定义E1/E2/nu1/nu2 % ========== 工况与求解参数 ========== Fr = 5000; % 径向载荷，N Fa = 3000; % 轴向载荷，N err_tol = 1e-6; % 位移收敛容差，归一化后无量纲 max_iter = 50; % NR最大迭代次数

单位约定是整个代码的命门：长度一律mm，力一律N，弹性模量一律N/mm²。这样K的单位自然就是N/mm^1.5，接触变形δ在mm量级（深沟球轴承通常0.01~0.1 mm），残差在10的负几次方，数值尺度比较均匀。如果你习惯用国际单位制，长度全换米，那K会变得非常大，迭代时数值梯度容易出现奇异——不是不能算，但没必要给自己找麻烦。

3.3 钢球载荷分布主循环

求解器主迭代的骨架如下，配合注释看逻辑：

% NR迭代求解五自由度位移向量 q = [x; y; z; thx; thy] q = zeros(5,1); % 初值：全部从零开始 load_steps = 10; % 载荷逐步施加的步数 for step = 1:load_steps F_ext_step = (step / load_steps) * F_ext; % 按比例加载 for iter = 1:max_iter [res, Qj, alpha_j] = compute_residual(q, F_ext_step, params); if norm(res) < err_tol * norm(F_ext_step) break; % 残差相对外部载荷足够小 end J = numerical_jacobian(@(q_) compute_residual(q_, F_ext_step, params), q); dq = -J \ res; % 标准NR步 q = q + dq; % 更新位移 end end

逻辑说明：

• compute_residual函数返回三个东西：五维残差向量res、钢球载荷数组Qj、接触角数组alpha_j。让Qj和alpha_j跟着走，是因为后处理画载荷分布图时还要用，避免重复计算。
• 外层load_steps是载荷逐步施加。刚开始迭代不建议直接把5000N全压上去，尤其是游隙较大的轴承，初始猜测q=0离真实解太远，NR容易直接发散。分成10步的话，第一步只加500N，解出来一个中间位移，作为下一步的初值——这也是工程里最常见的收敛处理手段。
• numerical_jacobian是数值差分求雅可比。核心代码在下一章展开，这里先理解整体闭环：给定位移算残差，用雅可比修正位移，残差够小就收敛。

跑完之后Qj数组就是Z个钢球的载荷分布，画出来应该是一个类似余弦钟形的曲线，最大载荷出现在最下方（径向载荷正下方）的那个钢球附近。如果画出来是锯齿形或不对称的形状，优先怀疑方位角起始点ψ0和坐标轴方向是否匹配（后续避坑章节会细说）。

4. Newton-Raphson迭代的工程设置：收敛判据、步长与初值

4.1 雅可比矩阵的数值差分技巧

准动力学模型的雅可比矩阵没有解析表达式，因为每个钢球的接触角变化对载荷的偏导写起来非常复杂。工程上直接数值差分，代价是每次迭代多算6次残差，但对5自由度模型来说开销完全可以接受。

function J = numerical_jacobian(fun, q) n = length(q); J = zeros(n); eps_step = 1e-7; % 差分步长，除以位移量级做归一化 scale = max(abs(q), 1e-3); % 防止某自由度为0时步长为0 for k = 1:n q_plus = q; q_minus = q; h = eps_step * scale(k); q_plus(k) = q_plus(k) + h; q_minus(k) = q_minus(k) - h; r_plus = fun(q_plus); r_minus = fun(q_minus); J(:,k) = (r_plus - r_minus) / (2 * h); % 中心差分 end end

关键在h = eps_step * scale(k)这行。如果直接用固定1e-7的绝对步长，而某个自由度的位移量级是1e-5，那差分步长比位移本身还大，雅可比列就是纯噪声。中心差分比单侧差分精度高一个数量级，代价是每列多算一次残差。5自由度模型每算一次雅可比需要10次残差调用，单次残差计算里包含Z个钢球的逐个接触判断，对Z=12的轴承整个耗时在0.1秒量级，完全不需要追求解析雅可比。

4.2 收敛判据与载荷逐步施加

收敛判据只用norm(res) < tol是不够稳的。位移的收敛和残差的收敛要同时看，尤其是载荷很大时，残差就算收敛到1e-4，位移可能还在缓慢漂移，对刚度矩阵输出影响明显。

我习惯用双判据：位移增量范数小于1e-7且残差相对值小于1e-6。代码块里加一条位移增量判断：

dq_norm = norm(dq); res_norm = norm(res) / norm(F_ext_step); if res_norm < err_tol && dq_norm < 1e-7 * (1 + norm(q)) break; end

注意dq_norm的判据要带(1 + norm(q))的缩放项——如果当前位移量级是0.1 mm，1e-7的绝对增量精度已经很高；但如果位移接近0（比如纯轴向受载时径向位移理论上是0），绝对判据反而容易提前终止。带上相对缩放项后，对零解附近的工况也稳。

载荷逐步施加的步数建议在5~10之间。并不是越多越好——load_steps设成100的话，每步载荷增量太小，前几步解出来的位移都接近0，白白浪费计算时间；设成3以下又起不到引导收敛的作用。实际调参时看迭代日志：如果前几步NR迭代次数都超过10次才收敛，说明步长太大；如果每步1~2次就收敛且残差立刻掉到容差内，说明步长太保守，可以减步数。

4.3 参数推荐速查表

整套参数下，常见的深沟球轴承和角接触球轴承都能在几十次迭代内收敛。如果你发现迭代卡在某个残差量级降不下去，不要急着调容差，先检查是不是雅可比计算里出现了某个钢球刚好在接触/脱离边界的位置——这种边界切换导致的残差不连续，是数值振荡的主要来源，下一章详细说。

5. 避坑与常见问题排查：五条血泪经验

5.1 接触刚度K算出来偏大三个数量级

现象：载荷分布图上每个钢球的Qj都异常大，最大载荷远超外部载荷总和，NR迭代一直不收敛或者收敛到明显错误的解。

原因：单位陷阱。查了参数表发现弹性模量直接用了GPa数值（比如207），没有换算成N/mm²（应该是207000）。K值因此差了约1000倍，接触变形被严重低估，所有钢球都显得刚性极大。

解决：在params_bearing.m开头加一行强制校验：

assert(E > 10000, '弹性模量E单位应为N/mm^2，当前值疑似使用了GPa');

从那以后我每次新建轴承参数都先跑一次断言，确认E在10万量级、钢球直径在毫米量级、载荷在千牛量级，三个数值尺度都对得上再进求解器。

5.2 载荷分布图不对称，峰值偏移到侧面

现象：轴向载荷为0、纯径向载荷作用下，载荷分布应该以载荷方向为对称轴。但画出来的Qj曲线明显偏向左或偏向右，峰值不在正下方。

原因：方位角ψj的起始位置和坐标轴方向不匹配。如果径向载荷Fr定义为沿x轴负方向（向下），第一个钢球ψ0在0度，那最大载荷位置应该在π方向附近的钢球上。如果ψ0设置成π/2，整个分布就像被旋转了90度。

解决：先看一眼代码里的Fr方向定义，再对照balance角计算。最省事的做法：用一个已知工况验证——纯径向载荷下的最大载荷对称性。如果不对称，就调整ψ0为0或π，直到左右钢球载荷关于ψ=π对称。这个检查我每次换轴承型号都会跑一遍，因为不同代码模板的坐标习惯不一样。

5.3 迭代残差卡在1e-2量级下不去

现象：NR迭代前半程收敛很快，残差从10降到0.1，然后卡在0.01~0.03之间振荡，继续迭代也不改善。

原因：某个钢球正好处于接触和脱离的临界状态。游隙大或轻载时，部分钢球的δj在0附近正负跳变，载荷从有到无的切换让残差函数出现折线拐点，导数在拐点处不连续，数值雅可比的梯度方向就乱了。

解决：给钢球载荷加一个光滑过渡，工程上叫正则化。当δj < 0.1 × 平均接触变形时，载荷不是直接置0，而是用二次曲线过渡到0：

if delta_j > 0 if delta_j < delta_cut Qj = K * delta_cut^1.5 * (delta_j / delta_cut)^2; % 光滑过渡 else Qj = K * delta_j^1.5; end else Qj = 0; end

delta_cut取平均接触变形的一小部分（比如1e-3 mm），物理含义是接触边缘的渐入渐出。加了这层过渡后，残差曲线不再有硬拐点，NR迭代就能稳定收敛到1e-8量级。如果还想保留“钢球要么接触要么脱开”的严格物理模型，那就接受残差只能到1e-3量级，但这会让后续刚度矩阵输出带噪声，我建议还是用光滑过渡。

5.4 高速工况下结果明显失真

现象：转速超过轴承极限转速的50%时，算出来的载荷分布和重载荷工况混在一起，区分不出离心力拉大内圈接触、拉小外圈接触的典型特征。

原因：准动力学模型默认忽略钢球离心力和陀螺力矩。低速时这两个量可以忽略，但高速时钢球离心力可达自身重力的几十倍，会显著改变内外圈接触载荷比——这是准动力学模型的适用边界。

解决：在compute_residual函数里加入离心力修正项。钢球公转转速nc = 0.5 × ni（外圈固定时），离心力为m × ωc² × dm/2。把它加到外圈接触载荷上，内圈载荷相应减小。修正后模型在中高转速下仍然可用，但如果转速超过DN值（节圆直径mm × 转速rpm）约50万，摩擦发热和弹流润滑效应开始主导，准动力学模型就不够用了，需要上全动力学或热-力耦合模型。代码注释里应该标明这个适用上限，我用下来觉得深沟球轴承在DN值30万以下这个模型和Romax的结果对得很准。

5.5 空载或反向受载时出现NaN

现象：外部载荷全部为0，或轴向载荷方向设反，程序报错或者位移解全是NaN。

原因：轴向位移和径向位移都把对应的外部载荷设成了0，但代码里某处对δj开根号——当δj = 0时0的1.5次方没问题，但当δj因为几何计算误差变成负数时，K × δj^1.5就是对负数开方，结果直接NaN。另外一个原因是纯轴向加载时径向方向的雅可比矩阵退化，因为径向刚度很小但也不是零，数值差分算出来可能接近奇异。

解决：在所有计算Qj的地方，先加delta_j > 0的判断再开方，这是最基本的健壮性防护。另外对雅可比矩阵求逆前，加一个小的对角正则化：

J_reg = J + 1e-6 * eye(5) * (1 + max(abs(J(:)))); dq = -J_reg \ res;

这个1e-6的相对正则化不会改变解的物理意义，但能防止矩阵奇异导致求解崩溃。从那以后我每次跑新工况都会先跑零载荷和反向载荷的边界测试，确认程序不会爆NaN再上正式工况。

6. 结果验证与后处理技巧：让模型可复现、可解释

跑通代码只是第一步，验证结果对得上物理常识才是真正把模型握在手里。最常用的检查是与经典Harris理论解的对照：最大钢球载荷Qmax ≈ 5Fr / (Z × cosα)（深沟球轴承α=0时简化为5Fr/Z），这个近似公式在纯径向载荷、小游隙条件下误差在10%以内。如果代码算出的Qmax偏离这个值超过15%，先查接触刚度K和游隙参数，不要怀疑公式——公式在低游隙下精度可靠。

后处理里我最常用的验证图有两个。第一个是极坐标载荷分布图，把每个钢球的Qj按方位角画成径向花瓣图，形状应该是光滑的单峰或双峰曲线，双峰出现在同时承受径向和轴向载荷时：

psi_deg = (0:Z-1) * 360 / Z; polarplot(deg2rad(psi_deg), Qj, 'o-'); title('钢球载荷分布（N）');

第二个是径向载荷-位移曲线。保持轴向载荷不变，Fr从0扫到额定动载荷的80%，记录径向位移x，画出来应该是一条斜率逐渐变陡的曲线——因为接触刚度随载荷增大而增大，这是Hertz非线性的直接体现。如果这条曲线是直线，说明模型里某个地方被线性化了，多半是K设成了常数而不是随曲率变化。

最后验证刚度矩阵输出的合理性：把位移解回代到总载荷公式，合成载荷应该和外部输入载荷一致，误差在1e-6量级，这是最直接的残差验证，也是程序自洽性的最终证明。

我从那以后每次算完一组工况，都会强制把等效刚度画成随载荷变化的曲线，看趋势、看数值量级、看和近似公式的偏差，三步都通过才敢把结果交给下游转子动力学分析。这套习惯帮我挡掉了至少三次单位错误和一次坐标方向错误。希望这份代码和笔记里的经验能帮你少走几趟弯路，直接拿到可信的轴承模型结果。

本文还有配套的精品资源，点击获取