振动力学学习笔记: 单自由度系统(三) 单自由度粘性阻尼系统的自由振动

前情提要
绪论(一) 振动力学的基本概念
绪论(二) 振动力学的基本问题与基本方法
绪论(三) 简谐振动及其三角函数、矢量、复数表示法

理想元件(一) 弹性元件及其简化模型
理想元件(二) 阻尼元件及其简化模型
理想元件(三) 质量元件及其简化模型

单自由度系统(一) 单自由度线性系统的数学模型
单自由度系统(二) 单自由度无阻尼系统的自由振动

---

---

本节讨论对象：单自由度粘性阻尼系统的自由振动

• 单自由度 (one degree of freedom) 系统在振动过程中任意瞬时的几何位置只需要一个独立坐标就能描述
• 粘性阻尼 (viscous damping) = 线性阻尼 与速度成正比的阻尼
• 自由振动 (free vibration) 振动系统仅受到初始条件 (初始位移、初始速度) 的激励，在后续的振动过程中没有其他激励

工程实际中的振动系统或大或小总有一些阻尼，因此，相比于无阻尼模型，粘性阻尼模型更加接近实际情况。

---

单自由度粘性阻尼系统自由振动的微分方程

通过受力分析，可以获得单自由度粘性阻尼系统自由振动的微分方程 $$m\ddot{x}(t)+c\dot{x}(t)+kx(t)=0\text{\hspace{1em}(1)}$$ 其中

• $m$ 为质量块质量
• $c$ 为粘性阻尼系数，单位为 $\text{N}\cdot \text{s}/\text{m}$
• $k$ 为弹簧刚度系数，单位为 $\text{N}/\text{m}$
  

将式 $(1)$ 方程两边同时除以质量 $m$，则方程可表示为 $$\ddot{x}(t)+2\zeta {\omega }_{\text{n}}\dot{x}(t)+{\omega }_{\text{n}}^{2}x(t)=0\text{\hspace{1em}(2)}$$ 其中

• 系统无阻尼时的固有频率 (单位 $\text{rad}/\text{s}$) $${\omega }_{\text{n}}=\sqrt{\frac{k}{m}}\text{\hspace{1em}(3)}$$
• 粘性阻尼比 (viscous damping ratio) = 粘滞阻尼因子 (viscous damping factor) (无量纲)$$\zeta =\frac{c}{2m{\omega }_{\text{n}}}=\frac{c}{2\sqrt{mk}}\text{\hspace{1em}(4)}$$

---

单自由度粘性阻尼系统自由振动解的三种情况

式 $(2)$ 是一个典型的二阶常系数齐次微分方程，其特征方程为 $$r^2+2\zeta {\omega }_{\text{n}}r+{\omega }_{\text{n}}^2=0\text{\hspace{1em}(5)}$$ 其特征根为 $$r_{1,2}=\left(-\zeta \pm \sqrt{{\zeta }^2-1}\right){\omega }_{\text{n}}\text{\hspace{1em}(6)}$$ 显然，根据阻尼比 $\zeta$ 的不同，特征根有不同的特点，进而产生不同的振动解，对应不同的振动类型，需进行分类讨论

---

小阻尼

$\zeta <1$ 的情况为 小阻尼 (underdamping)。此时特征根为一对共轭复数 $$r_{1,2}=\left(-\zeta \pm i\sqrt{1-{\zeta }^2}\right){\omega }_{\text{n}}=-\zeta {\omega }_n\pm i{\omega }_{\text{d}}$$ 其中 $${\omega }_{\text{d}}=\sqrt{1-{\zeta }^2}{\omega }_{\text{n}}\text{\hspace{1em}(7)}$$ 系统通解为 $$\begin{array}{rl}x(t)& =C_1{e}^{r_{1}t}+C_2{e}^{r_{2}t}=C_1{e}^{\left(-\zeta {\omega }_{\text{n}}+i{\omega }_{\text{d}}\right)t}+C_2{e}^{\left(-\zeta {\omega }_{\text{n}}-i{\omega }_{\text{d}}\right)t}\\ & =e^{-\zeta {\omega }_{\text{n}}t}\left[C_1\left(\cos {\omega }_{\text{d}}t+i\sin {\omega }_{\text{d}}t\right)+C_2\left(\cos {\omega }_{\text{d}}t-i\sin {\omega }_{\text{d}}t\right)\right]\\ & =e^{-\zeta {\omega }_{\text{n}}t}\left[\left(C_1+C_2\right)\cos {\omega }_{\text{d}}t+i\left(C_1-C_2\right)\sin {\omega }_{\text{d}}t\right]\end{array}$$ 可以看出，系统在平衡位置附近作往复振动，但振幅不断衰减。随 $C_1$ 和 $C_2$ 的变化，系统振动频率始终等于 ${\omega }_{\text{d}}$。因此，${\omega }_{\text{d}}$ 就是小阻尼系统的固有频率，称为 有阻尼固有频率 (damped natural frequency)

若令 $X_1=C_1+C_2$，$X_2=i\left(C_1-C_2\right)$，则小阻尼系统的通解可以表示为 $$x(t)=e^{-\zeta {\omega }_{\text{n}}t}\left(X_1\cos {\omega }_{\text{d}}t+X_2\sin {\omega }_{\text{d}}t\right)\text{\hspace{1em}(8)}$$ 式中 $X_1$ 和 $X_2$ 由初始条件确定。已知 $t=0$ 时系统的初位移为 $x_0$，初速度为 $v_0$，则有 $$\begin{array}{rl}x(t=0)& =e^{-\zeta {\omega }_{\text{n}}t}\left(X_1\cos {\omega }_{\text{d}}t+X_2\sin {\omega }_{\text{d}}t\right)=X_1=x_0\\ \dot{x}(t=0)& =e^{-\zeta {\omega }_{\text{n}}t}\left[\left(-\zeta {\omega }_{\text{n}}{X}_1+{\omega }_{\text{d}}{X}_2\right)\cos {\omega }_{\text{d}}t+\left(-\zeta {\omega }_{\text{n}}{X}_2-{\omega }_{\text{d}}{X}_1\right)\sin {\omega }_{\text{d}}t\right]\\ & =-\zeta {\omega }_{\text{n}}{X}_1+{\omega }_{\text{d}}{X}_2=v_0\end{array}$$ 由此可确定系数 $$\begin{array}{rl}{X}_1& =x_0\\ X_2& =\frac{v_0+\zeta {\omega }_{\text{n}}{x}_0}{{\omega }_{\text{d}}}=\frac{v_0+\zeta {\omega }_{\text{n}}{x}_0}{\sqrt{1-{\zeta }^2}{\omega }_{\text{n}}}\end{array}\text{\hspace{1em}(9)}$$ 若以 $X_1$ 和 $X_2$ 为直角边作直角三角形，记斜边为 $X$，$X_1$ 的对角为 ${\phi }_0$，则式 $(8)$ 可以被表示为 $$\begin{array}{rl}x(t)& =e^{-\zeta {\omega }_{\text{n}}t}\left(X_1\cos {\omega }_{\text{d}}t+X_2\sin {\omega }_{\text{d}}t\right)\\ & =e^{-\zeta {\omega }_{\text{n}}t}\left(X\sin {\phi }_0\cos {\omega }_{\text{d}}t+X\cos {\phi }_0\sin {\omega }_{\text{d}}t\right)\\ & =e^{-\zeta {\omega }_{\text{n}}t}\cdot X\sin \left({\omega }_{\text{d}}t+{\phi }_0\right)\end{array}\text{\hspace{1em}(10)}$$

如图为一小阻尼系统自由振动时的位移时程曲线，该系统质量 $m=1\text{kg}$，弹簧刚度系数 $k=1\text{N/m}$，粘性阻尼系数 $c=0.1\text{N}\cdot \text{s/m}$ (对应阻尼比 $\zeta =0.05$)，初始位移 $x_0=1\text{m}$，初始速度 $v_0=5\text{m/s}$

• 自由振动频率 = 有阻尼固有频率 $${\omega }_{\text{d}}=\sqrt{1-{\zeta }^2}{\omega }_{\text{n}}=\sqrt{\frac{\left(1-{\zeta }^2\right)k}{m}}$$ 可以看出，随着阻尼比的增大，系统的固有频率将会减小，阻尼比接近 $1$ 的时候，固有频率接近于 $0$
  

• 有阻尼固有周期 (damped natural cycle) $$T_{\text{d}}=\frac{2\pi }{{\omega }_{\text{d}}}=\frac{2\pi }{{\omega }_{\text{n}}\sqrt{1-{\zeta }^2}}=\frac{T_{\text{n}}}{\sqrt{1-{\zeta }^2}}\text{\hspace{1em}(11)}$$ 显然，对于给定系统，这是个定值。意味着阻尼的出现不会破坏系统自由振动的等时性 ，但因振幅有衰减，有阻尼系统不再具有周期性。其固有周期相比于无阻尼系统有所延长，阻尼越大延长越多。阻尼比接近 $1$ 时，有阻尼固有周期趋近于无穷大，这意味着下一个周期的开始大约在一 ~ 千年 ~ 以后——

• 振幅 $e^{-\zeta {\omega }_{\text{n}}t}\cdot X$ 按照指数规律衰减，其中 $$X=\sqrt{X_1^2+X_2^2}=\sqrt{x_0^2+\frac{{\left(v_0+\zeta {\omega }_{\text{n}}{x}_0\right)}^2}{\left(1-{\zeta }^2\right){\omega }_{\text{n}}^2}}\text{\hspace{1em}(12)}$$

• 对数衰减率 (logarithmic decrement) 一个自然周期相邻两个振幅之比的自然对数，反映振幅衰减的快慢 $$\delta =\ln \frac{e^{-\zeta {\omega }_{\text{n}}t}}{e^{-\zeta {\omega }_{\text{n}}(t+T_{\text{d}})}}=\zeta {\omega }_{\text{n}}{T}_{\text{d}}=\frac{2\pi \zeta }{\sqrt{1-{\zeta }^2}}\text{\hspace{1em}(13)}$$ 可以看出，振幅衰减的速度和阻尼大小呈现正相关。特别地，当阻尼很小，$1-{\zeta }^2\approx 1$，对数衰减率可以近似为 $$\delta \approx 2\pi \zeta \text{\hspace{1em}(14)}$$ 对数衰减率近似值相对于精确解的误差为 $$\epsilon =\frac{\left|2\pi \zeta -\frac{2\pi \zeta }{\sqrt{1-{\zeta }^2}}\right|}{\frac{2\pi \zeta }{\sqrt{1-{\zeta }^2}}}\times 100\%=\left(1-\sqrt{1-{\zeta }^2}\right)\times 100\%\text{\hspace{1em}(15)}$$ 当阻尼比小于 $0.31$ 时，误差小于 $5\%$，一般可以满足工程需要
  

• 初相位 (initial phase angle) ${\phi }_0$ $t=0$ 时的相位角，根据式 $(9)$ 有 $${\phi }_0=\arctan \frac{X_1}{X_2}=\arctan \frac{\sqrt{1-{\zeta }^2}{\omega }_{\text{n}}{x}_0}{v_0+\zeta {\omega }_{\text{n}}{x}_0}\text{\hspace{1em}(16)}$$

---

临界阻尼

刚刚在小阻尼的情况中，我们讲到，当阻尼比 $\zeta$ 接近 $1$ 时，振动频率趋近于 $0$，周期趋近于无穷大。而当阻尼比真的达到 $1$ ，即 $\zeta =1$ 时，就形成了 临界阻尼 (critical damping) 的情况。它介于小阻尼和后面要介绍的过阻尼之间

此时有 临界阻尼系数 $$c_0=2m{\omega }_n=2\sqrt{mk}\text{\hspace{1em}(17)}$$ 按照临界阻尼系数的概念，阻尼比也可以定义为 $$\zeta =\frac{c}{c_0}\text{\hspace{1em}(18)}$$

临界阻尼情况下，特征根为两个相同的负实数，即 $$r_1=r_2=-{\omega }_{\text{n}}$$ 系统通解为 $$x(t)=\left(C_1+C_{2}t\right)e^{-{\omega }_{\text{n}}t}\text{\hspace{1em}(19)}$$ 此时系统的振动按照指数规律衰减，不会产生周期震荡。式中 $C_1$ 和 $C_2$ 由初始条件确定。已知 $t=0$ 时系统的初位移为 $x_0$，初速度为 $v_0$，则有 $$\begin{array}{rl}x(t=0)& =\left(C_1+C_{2}t\right)e^{-{\omega }_{\text{n}}t}=C_1=x_0\\ \dot{x}(t=0)& =\left[-{\omega }_{\text{n}}\left(C_1+C_{2}t\right)+C_2\right]e^{-{\omega }_{\text{n}}t}=-{\omega }_{\text{n}}{C}_1+C_2=v_0\end{array}$$ 由此可确定系数 $$\begin{array}{rl}{C}_1& =x_0\\ C_2& =x_0{\omega }_{\text{n}}+v_0\end{array}\text{\hspace{1em}(20)}$$

如图为一临界阻尼系统自由振动时的位移时程曲线，该系统质量 $m=1\text{kg}$，弹簧刚度系数 $k=1\text{N/m}$，粘性阻尼系数 $c=2\text{N}\cdot \text{s/m}$ (对应阻尼比 $\zeta =1$)，初始位移 $x_0=1\text{m}$，初始速度 $v_0=5\text{m/s}$

---

过阻尼

$\zeta >1$ 的情况为 过阻尼 (overdamping)。此时特征根为两个不同的负实数，其中 $$\begin{array}{rl}{r}_1& =\left(-\zeta +\sqrt{{\zeta }^2-1}\right){\omega }_{\text{n}}\in (-{\omega }_n,0)\\ r_2& =\left(-\zeta -\sqrt{{\zeta }^2-1}\right){\omega }_{\text{n}}<-{\omega }_n\end{array}$$ 如下图，随阻尼比的增大，$r_1$ 越来越接近 $0$，$r_2$ 则越来越小

此时系统的通解为 $$x(t)=C_1{e}^{r_{1}t}+C_2{e}^{r_{2}t}\text{\hspace{1em}(21)}$$ 由于特征根为负，通解两项均按指数规律衰减。这是因为阻尼足够大，由初始激励输入给系统的能量很快就会被消耗掉，系统很快就会趋于平衡位置，而来不及产生反复的震荡。

$C_1$ 和 $C_2$ 均由初始条件，即初始时刻系统的位移和速度确定。已知 $t=0$ 时系统的初位移为 $x_0$，初速度为 $v_0$，则有 $$\begin{array}{rl}x(t=0)& =C_1+C_2=x_0\\ \dot{x}(t=0)& =C_1{r}_1{e}^{r_{1}t}+C_2{r}_2{e}^{r_{2}t}=C_1{r}_1+C_2{r}_2=v_0\end{array}$$ 由此可确定系数
$$C_{1,2}=\frac{x_0}{2}\left(1\pm \frac{\zeta {\omega }_{\text{n}}+v_0}{\sqrt{{\zeta }^2-1}{\omega }_{\text{n}}}\right)\text{\hspace{1em}(22)}$$

如图为一过阻尼系统自由振动时的位移时程曲线，该系统质量 $m=1\text{kg}$，弹簧刚度系数 $k=1\text{N/m}$，粘性阻尼系数 $c=5\text{N}\cdot \text{s/m}$ (对应阻尼比 $\zeta =2.5$)，初始位移 $x_0=1\text{m}$，初始速度 $v_0=5\text{m/s}$

---

生成单自由度线性系统自由振动时程曲线的MATLAB函数

根据上面计算的结果，编写如下函数。输入系统质量、弹性系数、粘性阻尼系数、初始条件、结束时间、步长，即可输出位移随时间的变化情况。如选择输出图像，还可直接生成位移时程曲线。

function [x,t] = x_free_vib(m,k,c,x0,v0,te,s,p)

%% Generating Displacement-time history data of linear vibration system with 1 DoF

%%   [OUTPUT] System charasteristics
%     x     Displacement data
%     t     Time data
%     a     Amplitude data

%%   [INPUT] System charasteristics
%     m     Mass
%     k     Elastic coefficient
%     c     Damping coefficient

%%   [INPUT] Initial conditions 
%     x0    Initial displacement
%     v0    Initial velocity

%%   [INPUT] Calculating settings 
%     te    End time
%     s     Steps (Default=100)
%     p     Plot or not (Default=0=No)


if nargin<7
    s=100;% Set default step
end
if nargin<8
    p=0;
end

t=0:te/s:te;% Generating time series
wn=sqrt(k/m); % Calculating undamped natural frequency
z=c/(2*sqrt(m*k));% Calculating damping ratio to determine formula of the general solution
if z<1 % Underdamping
    wd=sqrt(1-z^2)*wn;% Calculating damped natural frequency
    X1=x0;
    X2=(v0+z*wn*x0)/(wd);
    x=exp(-z*wn*t).*(X1*cos(wd*t)+X2*sin(wd*t));
    a=exp(-z*wn*t)*sqrt(X1^2+X2^2);
elseif z==1 % Critical damping
    C1=x0;
    C2=x0*wn+v0;
    x=(C1+C2*t).*exp(-wn*t);
else % Overdamping
    r1=(-z+sqrt(z^2-1))*wn;
    r2=(-z-sqrt(z^2-1))*wn;
    C1=x0/2*(1+(z*wn+v0)/(sqrt(z^2-1)*wn));
    C2=x0/2*(1-(z*wn+v0)/(sqrt(z^2-1)*wn));
    x=C1*exp(r1*t)+C2*exp(r2*t);
end
if p
    plot(t,x)
        axis([0,max(t),-ceil(max(abs(x))*1.1),ceil(max(abs(x))*1.1)])
    if z<1
        plot(t,a,t,-a,'LineStyle','--')
    end
end
end

调用语法

Matlab 代码	注释
x = x_free_vib(m,k,c,x0,v0,te)	返回系统位移值的时间序列 m 质量 | k 弹性系数 | c 阻尼系数 x0 初始位移 | v0 初始速度 te 结束时间 | 默认分析步数为100
x = x_free_vib(m,k,c,x0,v0,te,s)	还指定分析步数 s
[x,t] = x_free_vib(m,k,c,x0,v0,te,s)	还返回时间的序列 t
[x,t] = x_free_vib(m,k,c,x0,v0,te,s,1)	还绘制图像
x_free_vib(m,k,c,x0,v0,te,s,1);	只绘制图像，不返回数据

---

不同阻尼系数下的系统自由振动的位移时程曲线

使用上一节编写的函数，绘制如下六个系统的位移时程曲线，探究阻尼系数对系统振动特性的影响

系统编号		1	2	3	4	5	6
		(无阻尼)	(小阻尼)		(临界阻尼)	(过阻尼)
系统特性	质量	1
	弹簧刚度系数	1
	粘性阻尼系数	0	0.2	1	2	3	5
	阻尼比	0	0.1	0.5	1	1.5	2.5
初始条件	初始位移	1
	初始速度	5
计算设置	终止时间	60
	步数	1000

MATLAB 程序

clear all
figure
hold on
x_free_vib(1,1,5,1,5,60,1000,1);
x_free_vib(1,1,3,1,5,60,1000,1);
x_free_vib(1,1,2,1,5,60,1000,1);
x_free_vib(1,1,1,1,5,60,1000,1);
x_free_vib(1,1,0.2,1,5,60,1000,1);
x_free_vib(1,1,0,1,5,60,1000,1);

运行结果

注：为使图片更加清晰美观，使用图像编辑软件对线型、线宽、标记等进行了一些处理。这些处理也可在 MATLAB 中完成，用哪个取决于个人喜好。非重点，不赘述。

从图中可以看出，对于 $\zeta <1$ 的小阻尼系统，随着阻尼的增大，振幅衰减的速度提升，同时震荡周期发生一定程度的减小。

阻尼比增加到 $1$，系统进入临界阻尼状态，震荡现象消失。

阻尼比进一步增大，$\zeta >1$ 时，系统成为过阻尼系统。随着阻尼的增大，系统位移的峰值减小。值得注意的是，并非是阻尼越大衰减越快。图中衰减最快的是临界阻尼的情况。

实际工程中，我们有时关注的是系统位移的峰值，有时关注的是衰减到给定值的速度，更多时候两者都要关注；此外，还有能耗、一致性等多方面的因素需要兼顾。因此，阻尼比不是越大越好，也不是越小越好，需要我们根据实际情况进行选择。

---

参考文献

[1] 鲍文博,白泉,陆海燕.振动力学基础与MATLAB应用[M].北京:清华大学出版社,2015:47~58.