1 引言

你是否注意过，提琴弓子摩擦琴弦能发出悠扬乐音，而推门时轴承却会发出刺耳的“吱嘎”声？这两种现象背后，其实隐藏着同一个物理机理——干摩擦自激振动。
在机械工程中，这类自振常表现为进给机构低速爬行、机床车刀切削颤振、制动器啸叫等有害现象，会大幅降低加工精度、加剧零部件磨损，严重时直接缩短设备服役寿命。

2 物理模型与数学描述

2.1 系统构成

选取经典滑块-弹簧-匀速传送带单自由度动力学模型：
质量为 $m$ 的滑块通过刚度 $k$ 的弹簧连接固定端，滑块放置在以恒定速度 $v_0$ 匀速运动的传送带上；滑块与传送带间存在非线性干摩擦力 $\varphi(v)$，摩擦力由滑块与传送带的相对速度决定。
定义 $\xi$ 为弹簧伸长量，则滑块与传送带的相对速度：

在这里插入图片描述

2.2 干摩擦力非线性特性

干摩擦力随相对速度变化具备典型Stribeck特征：

  1. 相对速度为0时，摩擦力达到最大静摩擦力 $F_s$；
  2. 滑块发生相对滑动瞬间，摩擦力骤降为库仑滑动摩擦力 $F_c$；
  3. 随着相对速度继续增大，摩擦力受粘性阻尼作用缓慢上升。

零速度附近摩擦力负斜率特性是干摩擦系统产生自激振动的核心能量来源，非线性阻尼会持续向系统输入振动能量。
在这里插入图片描述

2.3 系统动力学运动方程

为简化理论分析，对系统参数做归一化处理：$m=1，k=1$。
滑块受弹簧恢复力 $-\xi$、非线性干摩擦力 $\varphi(\dot{\xi}-v_0)$ 作用，由牛顿第二定律可得系统原始运动微分方程：

系统静平衡位置 $\xi_s$ 满足受力平衡条件：

引入坐标平移变换消除静平衡偏移：令 $x = \xi - \xi_s$，将平衡位置平移至坐标原点，定义等效非线性阻尼函数：

代入原方程得到以平衡位置为原点的标准振动方程：

阻尼函数特性：在 $y=0$（低速区间）附近 $\Phi(y)$ 斜率为负，系统表现为负阻尼，振动能量持续累积；当速度绝对值较大时，阻尼斜率由负转正，系统正向耗散振动能量。

3 相平面分析与稳定极限环机理

3.1 一阶状态相空间方程

令状态变量 $x_1=x，x_2=\dot{x}=y$，将二阶微分方程改写为一阶自治微分方程组：

相轨迹微分形式：

3.2 零斜率等倾线与平衡点稳定性

在这里插入图片描述
零斜率等倾线（$dy/dx=0$）满足：

  1. 原点附近等倾线分布在一、三象限，平衡点为不稳定焦点；微小扰动下相轨迹螺旋向外发散，振动幅值不断增大，负阻尼持续向系统注入能量；
  2. 当振动幅值提升至高速区间，系统切换为正阻尼，振动能量开始耗散；
  3. 能量注入速率与能量耗散速率动态平衡时，相轨迹不再发散或收敛，最终闭合形成稳定极限环。

3.3 极限环物理含义：粘-滑周期性振荡

极限环对应机械系统典型粘滞-滑动（Stick-Slip）周期性运动：

  1. 粘滞阶段：滑块与传送带无相对滑动，弹簧持续拉伸积蓄弹性势能；
  2. 滑动触发：弹簧弹力超过最大静摩擦力，滑块相对传送带发生滑移；
  3. 滑移减速阶段：库仑摩擦力作用下滑块速度降低，弹簧释放能量；
  4. 再次粘滞：滑块速度衰减至与传送带同步，再次进入粘滞状态，循环往复形成等幅自激振荡。

4 MATLAB数值仿真实现

4.1 Stribeck干摩擦数学模型

采用工程通用Stribeck摩擦模型表征干摩擦非线性特性，为规避符号函数 $\text{sgn}(v)$ 在零速度处不连续导致的数值奇异，使用 $\tanh(\alpha v)$ 平滑近似符号函数：

仿真参数配置表

参数 物理含义 取值
$F_c$ 库仑滑动摩擦力 0.4
$F_s$ 最大静摩擦力 0.8
$v_s$ Stribeck特征速度 0.2
$\sigma$ 粘性摩擦系数 0.05
$\alpha$ 零速度平滑系数 100
$v_0$ 传送带匀速运行速度 0.8

4.2 完整MATLAB仿真代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
%% 单自由度干摩擦自振系统极限环 (教材模型)
clear; clc; close all;

% ---------- 参数 ----------
m = 1;          % 质量
k = 1;          % 弹簧刚度
v0 = 0.8;       % 平台速度

% 摩擦力参数 (Stribeck)
Fc = 0.4;       % 库仑摩擦
Fs = 0.8;       % 静摩擦
vs = 0.2;       % Stribeck速度
sigma = 0.05;   % 粘性系数
alpha = 100;     % tanh平滑因子

% 摩擦力函数 phi(v)
phi = @(v) Fc*tanh(alpha*v) + (Fs-Fc)*exp(-(v/vs).^2).*tanh(alpha*v) + sigma*v;

% 定义 Phi(y) = phi(y - v0) - phi(-v0)
Phi = @(y) phi(y - v0) - phi(-v0);

% ---------- 状态方程 (x1 = x, x2 = y = dx/dt) ----------
% 方程: dx1/dt = x2, dx2/dt = -x1 - Phi(x2)
odefun = @(t, z) [z(2); -z(1) - Phi(z(2))];

% ---------- 仿真设置 ----------
Tsim = 100;         % 足够长，进入稳态
dt = 0.001;
tspan = 0:dt:Tsim;
z0 = [0.45; 0.45];    % 小扰动 (原点不稳定，将发散到极限环)

options = odeset('RelTol', 1e-6, 'AbsTol', 1e-8, 'MaxStep', dt);
[t, Z] = ode45(odefun, tspan, z0, options);   % 此处可用ode45，因为非刚性

% 提取变量
x = Z(:,1);
y = Z(:,2);    % 速度

% ---------- 绘图 ----------

% 图1：相平面 (x vs y)
figure('Name', '相平面极限环');
plot(x, y, 'b-', 'LineWidth', 1.5);
xlabel('位移 x');
ylabel('速度 y = \dot{x}');
title('干摩擦自振系统极限环');
grid on; axis equal; hold on;
plot(x(1), y(1), 'ro', 'MarkerSize', 8, 'LineWidth', 2);
legend('相轨迹', '起始点');
hold off;

% 图2：摩擦力曲线 Phi(y) 与零等倾线
figure('Name', '摩擦力特性');
% 绘制理论曲线 Phi(y)
y_plot = linspace(-2, 2, 300);
Phi_plot = Phi(y_plot);
plot(y_plot, Phi_plot, 'r-', 'LineWidth', 1.5);
xlabel('速度 y');
ylabel('\Phi(y)');
title('非线性阻尼函数 \Phi(y)');
grid on;
hold off;

% 图3：时间响应 (可选)
figure('Name', '时间响应');
subplot(2,1,1);
plot(t, x, 'b-');
xlabel('时间 (s)'); ylabel('位移 x');
title('位移响应'); grid on;

subplot(2,1,2);
plot(t, y, 'r-');
xlabel('时间 (s)'); ylabel('速度 y');
title('速度响应'); grid on;

4.3 仿真结果分析

  1. 相平面极限环结果：系统从微小初始扰动出发，相轨迹由内向外螺旋发散，最终收敛至闭合曲线，形成全局稳定极限环。
    在这里插入图片描述

  2. 非线性阻尼特性：$\Phi(y)$ 在零速度附近斜率为负，负阻尼持续激励系统振动；高速区间阻尼转正，限制振幅无限增大，能量动态平衡是极限环产生的根本原因。
    在这里插入图片描述

  3. 时域响应特性：位移、速度时序曲线经过短暂过渡过程后变为等幅周期振荡，无衰减、无发散，对应机械结构持续低频自激振动。
    在这里插入图片描述