简介：这份资源面向从事机械传动、伺服控制与摩擦补偿研究的工程师及高校学生，提供Stribeck摩擦力模型的Matlab实现代码与调用示例，用于解决低速工况下负斜率摩擦现象的建模与仿真问题。压缩包共3个文件，均为m脚本，体积约1KB，分别承担模型参数定义、摩擦曲线计算与演示调用等职责，结构紧凑、便于直接嵌入既有仿真流程。资源围绕指数型摩擦模型展开，即摩擦力在克服静摩擦后先随速度增大而减小、随后再上升，可帮助读者复现Stribeck曲线并观察不同参数对负斜率区间的影响。目前已有1137人学习下载，适合作为摩擦建模入门与参数调试的参考脚本，读者可据此快速搭建仿真验证环境，理解静摩擦、库仑摩擦与速度阻尼之间的耦合关系，并在此基础上扩展自适应补偿或辨识算法。
1. Stribeck 摩擦模型：为什么速度趋近于零时摩擦力反而最难算
如果你在 Matlab 里搭过一个带摩擦的机械系统仿真，大概率遇到过这种怪事：速度过零附近，求解器突然把步长压到极小，仿真时间几乎停滞，或者干脆报出刚性方程警告。很多人第一反应是求解器参数没调好，换 ode45 换 ode15s，折腾一圈发现该卡还是卡。问题的根子不在求解器，而在摩擦模型本身——尤其是 Stribeck 曲线在零速附近那段陡峭的负斜率。
干摩擦、粘性摩擦、Stribeck 效应叠在一起，构成了工程上最常用的摩擦力描述。它的核心价值在于：能同时刻画静摩擦、边界润滑、混合润滑到流体润滑的过渡，比单纯的库仑摩擦加粘性项更贴近真实导轨、轴承、气缸的 behavior。但代价是，这个模型在零速附近是强非线性、甚至不连续的，直接丢给定步长或显式求解器，翻车是迟早的事。这份 Matlab 实现代码和使用示例，解决的正是「模型怎么建、零速附近怎么处理、仿真怎么跑稳」这三个层层递进的问题。搞懂它，你再看任何摩擦补偿、伺服控制、机械臂关节建模的论文，都不会被那堆公式唬住。
2. Stribeck 模型的数学结构：四段摩擦状态是怎么拼出来的
2.1 从库仑摩擦到 Stribeck 曲线：多出来的那一块是什么
最原始的摩擦力模型只有一项：F = μ·N，方向与运动趋势相反。这个模型在速度较大时勉强能用，但它解释不了一个现象——为什么物体从静止开始推动的瞬间，需要的力往往比匀速滑动时更大。静摩擦系数大于动摩擦系数，这是第一层修正。
第二层修正来自润滑。当接触面之间有油膜或边界润滑剂时，速度从零开始增加，摩擦力并不是单调上升，而是先下降再上升。先下降那段，就是 Stribeck 效应：低速时润滑膜尚未完全建立，微凸体接触多，摩擦大；速度增加到某个值，流体动压效应起来，接触被油膜隔开，摩擦降到最低；速度再高，粘性剪切占主导，摩擦又随速度线性上升。
所以完整的 Stribeck 模型通常写成四项之和：
% Stribeck 摩擦力模型核心表达式 % F_friction = Fc + (Fs - Fc) * exp(-(v/vs)^delta) + Fv * v + Foffset % Fc: 库仑摩擦力 (N) % Fs: 最大静摩擦力 (N) % vs: Stribeck 速度 (m/s)，决定过渡段宽度 % delta: Stribeck 形状因子，通常取 0.5~2 % Fv: 粘性摩擦系数 (N·s/m) % Foffset: 偏置项，用于补偿传感器零漂或预紧力 function F = stribeck_friction(v, Fc, Fs, vs, delta, Fv, Foffset) % v 为相对速度，标量或向量 F = Fc + (Fs - Fc) .* exp(-abs(v./vs).^delta) + Fv .* v + Foffset; end这段代码里，Fc + (Fs - Fc)*exp(...)描述的是速度相关的那部分摩擦，它在 v=0 时等于 Fs，在 |v| 远大于 vs 时趋近于 Fc。Fv*v是粘性项，Foffset是常数偏置。注意这里用了abs(v)，因为 Stribeck 项是速度大小的函数，而粘性项保留符号，这样摩擦力方向始终与速度方向一致。
参数vs是最关键的：它决定了摩擦力从静摩擦降到库仑摩擦的过渡速度范围。vs 越小，过渡越陡，零速附近越难算；vs 越大，曲线越平滑，但可能偏离实际物理。delta控制过渡曲线的形状，delta=1 时是指数衰减，delta=2 时是高斯型衰减，工程上常取 1 或 2。
2.2 零速附近的数学病态：为什么显式求解器会崩
把上面的表达式画出来，你会看到在 v=0 附近，摩擦力对速度的导数非常大。对exp(-|v/vs|^delta)求导，当 delta=1 时，导数在 v=0 处是 - (Fs-Fc)/vs * sign(v)，这是一个有限但很大的值；当 delta<1 时，导数在 v=0 处趋于无穷。这意味着摩擦力的雅可比矩阵在零速附近条件数极差，显式求解器为了维持稳定，只能把步长压到极小。
更麻烦的是，真实的摩擦还有「粘滞」状态：当驱动力小于最大静摩擦力时，物体速度保持为零，摩擦力大小等于驱动力，方向相反。这是一个约束条件，不是微分方程。很多仿真代码在这里用「速度阈值」来切换：|v| < v_th 时认为处于粘滞状态，摩擦力按静摩擦处理；否则按 Stribeck 曲线算。这个切换本身引入了不连续，求解器在切换点附近会反复试探，进一步拖慢仿真。
常见的处理思路有三种：一是用连续函数近似粘滞状态，比如把摩擦力写成F = Fs * tanh(k*v)加上 Stribeck 项，用 tanh 的平滑过渡代替硬切换；二是用事件检测，在速度过零时让求解器停下来重新初始化；三是直接换用隐式求解器 ode15s 或 ode23t，让求解器自己处理刚性。这份 Matlab 代码里，我一般会看到它同时提供了连续近似和事件检测两种模式，让使用者根据精度和速度需求选。
2.3 参数辨识：vs 和 delta 怎么定
Stribeck 模型的参数不是拍脑袋来的。Fc、Fs 可以从实验的力-位移曲线里读：匀速段平均力是 Fc，启动瞬间峰值力是 Fs。Fv 可以从高速段的斜率拟合。难的是 vs 和 delta，它们藏在过渡段里。
一个实用的做法是：先固定 delta=1 或 2，然后用最小二乘拟合 vs。把实验数据里的速度 v 和摩擦力 F 拿出来，构造残差F_meas - Fc - Fv*v - Foffset，这个残差应该等于(Fs-Fc)*exp(-(v/vs)^delta)。两边取对数，变成线性问题，可以直接解出 vs。如果 delta 也要辨识，就网格搜索 delta，对每个 delta 做上面的线性拟合，取残差最小的组合。
代码里通常会有一个fit_stribeck_params.m或者类似的函数，输入实验数据矩阵，输出参数结构体。注意实验数据要覆盖低速段，否则 vs 辨识不准。如果实验台最低只能跑到 0.01 m/s，而实际 vs 是 0.001 m/s，那拟合出来的 vs 会偏大，仿真时零速附近的摩擦特性就失真了。
3. Matlab 实现细节：从方程到可运行代码的几层封装
3.1 状态空间建模：把摩擦力塞进 ODE 右边
Matlab 里做动力学仿真，标准套路是写一个返回状态导数的函数，交给 ode45 或 ode15s。对于单自由度系统 m*x'' = F_drive - F_friction(x')，状态向量取 [x; v]，导数就是 [v; (F_drive - F_friction(v))/m]。
% 单自由度系统带 Stribeck 摩擦的 ODE 函数 function dydt = sdof_stribeck(t, y, params) x = y(1); v = y(2); % 驱动力，这里用正弦激励举例 F_drive = params.A * sin(params.omega * t); % 摩擦力计算 F_fric = stribeck_friction(v, params.Fc, params.Fs, ... params.vs, params.delta, params.Fv, params.Foffset); % 状态导数 dydt = [v; (F_drive - F_fric) / params.m]; end这个函数本身很简单，关键在于stribeck_friction里对零速的处理。如果直接用原始公式，v=0 时abs(v./vs).^delta是 0，exp 是 1，摩擦力等于 Fs，没问题。但 v 从正变负经过零点的瞬间，粘性项Fv*v变号，摩擦力方向突变，导数不连续。ode45 会在这里把步长砍到 1e-12 级别，仿真时间稍微长一点就跑不动。
改进方式是在stribeck_friction里加一个平滑因子。比如把abs(v)换成sqrt(v^2 + eps_smooth^2)，eps_smooth 取一个很小的正数，比如 1e-4。这样 v=0 附近函数变成光滑的，导数连续，求解器不用反复试探。代价是零速附近的摩擦力略微偏离理论值，但工程上完全可以接受。
3.2 粘滞状态的处理：阈值切换 vs 连续近似
粘滞状态是摩擦仿真里最容易被忽略的坑。当驱动力小于 Fs 时，物体应该不动，速度恒为零。但 ODE 求解器不知道这个约束，它会根据(F_drive - F_fric)/m算加速度，如果 F_fric 取 Fs，加速度可能是负的，速度变成负的，下一步摩擦力又反向，来回振荡。
代码里常见的处理是加一个速度阈值 v_eps，比如 1e-6 m/s。当 |v| < v_eps 且 |F_drive| < Fs 时，强制令 dydt(2) = 0，并且摩擦力等于 F_drive。这相当于在 ODE 右边做了一个条件判断，虽然破坏了光滑性，但物理上是对的。
% 带粘滞判断的摩擦力计算 function F = stribeck_friction_with_stiction(v, F_drive, params) v_eps = 1e-6; % 粘滞速度阈值 if abs(v) < v_eps && abs(F_drive) < params.Fs % 粘滞状态：摩擦力与驱动力平衡 F = F_drive; else % 滑动状态：按 Stribeck 曲线计算 F = params.Fc + (params.Fs - params.Fc) * ... exp(-abs(v/params.vs)^params.delta) + ... params.Fv * v + params.Foffset; end end这个写法在 v_eps 附近仍然有不连续，但比不加判断要好得多。如果还是卡，可以把 v_eps 适当放大，比如 1e-4，或者用 tanh 平滑：F = Fs * tanh(v/v_eps) + ...，这样连判断都省了，全连续。
3.3 仿真脚本的组织：参数、求解、绘图分离
一份能用的 Matlab 代码，不会把所有东西塞在一个脚本里。常见的组织方式是：
params.m：定义所有物理参数和仿真参数，返回一个结构体stribeck_friction.m：摩擦力计算函数sdof_stribeck.m：ODE 函数run_simulation.m：主脚本，调用 ode15s，画图
主脚本里关键的一行是求解器选择：
% 主仿真脚本片段 params = params(); y0 = [0; 0]; % 初始位移和速度 tspan = [0 10]; % 用 ode15s 处理刚性，或者 ode45 加事件检测 options = odeset('RelTol', 1e-6, 'AbsTol', 1e-8, ... 'MaxStep', 1e-3); [t, y] = ode15s(@(t,y) sdof_stribeck(t, y, params), tspan, y0, options); % 提取速度和摩擦力 v = y(:,2); F_fric = arrayfun(@(vi) stribeck_friction(vi, params.Fc, params.Fs, ... params.vs, params.delta, params.Fv, params.Foffset), v); % 画相图 figure; plot(v, F_fric, 'b-', 'LineWidth', 1.5); xlabel('速度 (m/s)'); ylabel('摩擦力 (N)'); title('Stribeck 摩擦曲线'); grid on;MaxStep设成 1e-3 是为了防止求解器在平滑段跳太大步，错过零速附近的细节。RelTol和AbsTol比默认值严一些，保证低速段精度。如果仿真还是慢，可以把MaxStep放宽到 1e-2，同时把vs适当调大，牺牲一点精度换速度。
绘图部分，除了相图，通常还会画时域响应：位移-时间、速度-时间、摩擦力-时间。这三张图放在一起，能直观看出粘滞-滑动切换的时刻，以及 Stribeck 效应在启动阶段的表现。
4. 避坑与排查：Stribeck 摩擦仿真里最常见的五个翻车点
4.1 现象：仿真在零速附近卡死，步长被压到 1e-15
原因：摩擦力表达式在 v=0 处导数不连续，或者粘滞判断引入了硬切换，求解器为了满足误差容限反复试探。
解决：把abs(v)换成sqrt(v^2 + eps^2)，eps 取 1e-4 到 1e-6；或者把粘滞判断改成 tanh 平滑。如果还卡，换 ode15s 并放宽RelTol到 1e-4。
4.2 现象：速度在零附近高频振荡，正负来回跳
原因：粘滞状态没有正确处理，求解器在F_drive - F_fric变号时来回过冲。
解决：加粘滞判断，当 |v| < v_eps 且 |F_drive| < Fs 时强制加速度为零。v_eps 不要取太小，1e-6 到 1e-4 之间试。如果振荡还在，检查 Fs 是否小于 Fc，那会导致模型自相矛盾。
4.3 现象：辨识出的 vs 偏大，低速段摩擦特性失真
原因：实验数据没有覆盖足够低的速度段，拟合时 vs 被高速段数据带偏。
解决：实验时尽量跑到最低速，比如 0.001 m/s。如果设备限制，可以在拟合时给低速数据加权，或者固定 delta 为 1，只拟合 vs，减少自由度。
4.4 现象：仿真结果对 delta 极其敏感，delta 从 1 改到 1.1 曲线就大变
原因：delta 控制过渡段形状，在 vs 附近，delta 的微小变化会显著改变 exp 项的衰减速率。
解决：delta 不要连续调，按物理意义选：边界润滑主导取 1，混合润滑取 1.5 到 2。如果非要辨识，用网格搜索而不是梯度下降，避免陷入局部极小。
4.5 现象：ode45 跑得动，ode15s 反而报错
原因：ode15s 是隐式求解器，需要雅可比矩阵。如果摩擦力函数里有if判断或abs，雅可比不连续，ode15s 的数值雅可比会算错。
解决：用 ode15s 时，摩擦力函数必须连续可导。把粘滞判断改成 tanh 平滑，把abs(v)改成sqrt(v^2+eps^2)。或者给 ode15s 提供解析雅可比，但 Stribeck 模型的雅可比推导比较繁琐，一般不值得。
5. 从代码到设计取舍：Stribeck 模型什么时候该用，什么时候该换
Stribeck 模型不是万能的。它的优势在于物理意义清晰、参数少、能覆盖从静摩擦到流体润滑的过渡。但它的局限也很明显：只适用于稳态摩擦，不能描述预滑动位移（Dahl 模型、LuGre 模型才能）；参数辨识依赖实验数据，vs 和 delta 的准确性直接决定仿真可信度；零速附近的数值处理永远是在精度和稳定性之间做妥协。
我一般会这么判断：如果仿真关注的是宏观运动，比如气缸往复、导轨进给，速度大部分时间远离零，Stribeck 足够用，而且计算量小。如果关注的是微米级定位、爬行现象、预滑动阶段的力-位移关系，那就得上 LuGre 或 GMS 模型，Stribeck 只能作为其中的稳态摩擦项。
另一个取舍是连续近似 vs 事件检测。连续近似代码简单，求解器友好，但零速附近的摩擦力是「糊」的，粘滞状态被平滑掉了。事件检测能精确捕捉粘滞-滑动切换时刻，但代码复杂，每次事件触发都要重新初始化求解器，仿真速度慢。如果做控制算法验证，需要精确的切换时刻，事件检测值得；如果只是看整体响应趋势，连续近似就够了。
还有一个容易被忽略的点：温度。Stribeck 参数，尤其是 Fv 和 vs，对温度敏感。高温下润滑剂粘度下降，Fv 减小，vs 可能偏移。如果仿真工况温度变化大，要么把参数写成温度的函数，要么在多个温度点分别辨识参数，仿真时插值。这份代码里如果没有温度模块，用的时候要心里有数。
从那以后我每次拿到一个新的摩擦仿真任务，都强制走一遍：先画实验的力-速度曲线，确认 Stribeck 效应是否明显；再检查零速附近的数据密度，决定 vs 能不能辨识准；最后在代码里把平滑因子和粘滞阈值显式写出来，不依赖默认值。这三步走完，仿真跑崩的概率会低很多。希望帮到你。
本文还有配套的精品资源，点击获取
