弹簧-高质量-阻尼实验指导书

质量-弹簧-阻尼系统实验教学指导书

北京理工大学机械与车辆学院

2016.3

实验一：单自由度系统数学建模及仿真 1 实验目的

（1）熟悉单自由度质量-弹簧-阻尼系统并进行数学建模； （2）了解MATLAB 软件编程，学习编写系统的仿真代码； （3）进行单自由度系统的仿真动态响应分析。 2 实验原理

单自由度质量-弹簧-阻尼系统，如上图所示。由一个质量为m 的滑块、一个

刚度系数为k 的弹簧和一个阻尼系数为c 的阻尼器组成。系统输入：作用在滑块上的力f (t )。系统输出：滑块的位移x (t )。

建立力学平衡方程：

m x c x kx f ???

++=

变化为二阶系统标准形式：

22f

x x x m

ζωω??

?

++=

其中：ω是固有频率，ζ是阻尼比。

ω=

2c m ζω=

= 2.1 欠阻尼(ζ<1)情况下，输入f (t )和非零初始状态的响应：

()()sin()))]

t t x t t d e ζωττζωττ

+∞

--=?

-=

-+-?

2.2 欠阻尼(ζ<1)情况下，输入f(t)=f0*cos(ω0*t) 和非零初始状态的的响应：

022

3

00

22222

00

222222

2

()cos(arctan())

2f

[(0)]cos()

[()(2)]

sin(

t

t

x t t

x e

k

e

ζω

ζω

ζωω

ω

ωω

ζωω

ωωζωω

-

?

-

=-

-

++

-+

+)

输出振幅和输入振幅的比值：A=

3 动力学仿真

根据数学模型，使用龙格库塔方法ODE45求解，任意输入下响应结果。

仿真代码见附件

4 实验

4.1 固有频率和阻尼实验

（1）将实验台设置为单自由度质量-弹簧-阻尼系统。

（2）关闭电控箱开关。点击setup菜单，选择Control Algorithm，设置选择Continuous Time Control，Ts=0.0042，然后OK。

（3）点击Command菜单，选择Trajectory，选取step,进入set-up,选取Open Loop Step设置(0)counts, dwell time=3000ms,(1)rep, 然后OK。此步是为了使控制器得到一段时间的数据，并不会驱动电机运动。

（4）点击Data菜单，选择Data Acquisition,设置选取Encoder#1 ，然后OK离开；从Utility菜单中选择Zero Position使编码器归零。

（5）从Command菜单中选择Execute，用手将质量块1移动到2.5cm左右的位置（注意不要使质量块碰触移动限位开关），点击Run, 大约1秒后，放开手使其自由震荡，在数据上传后点击OK。

（6）点击Plotting菜单，选择Setup Plot，选取Encoder #1 Position；然后点击Plotting菜单，选择Plot Data，则将显示质量块1的自由振动响应曲线。（7）在得到的自由振动响应曲线图上，选择n个连续的振幅明显的振动周期，计算出这段振动的时间t，由n/t即可得到系统的频率，将Hz转化为rad/sec即为系统的振动频率ω。

（8）在自由振动响应曲线图上，测量步骤7选取时间段内初始振动周期的振幅X0以及末尾振动周期的振幅Xn。由对数衰减规律即可求得系统阻尼比。

（9）实验数据记录

（10）在仿真代码基础上，计算出实验结果对应的理论结果。对比分析理论和实验结果的差异。完成实验报告。

4.2 幅频特性实验

（1）点击Command菜单，选择Trajectory，选取Sinuscidal，进入set-up,选取Open LoopStep设置Amplitude(0.5V), Frequency(2Hz),Repetition(8)，然后OK。（2）从Utility菜单中选择Zero Position使编码器归零。从Command菜单中选择Execute，点击Run,在数据上传后点击OK。

（3）然后点击Plotting菜单，选择Plot Data，则将显示滑块的受迫振动响应曲线。在响应曲线图上，测量出振动振幅，计算出振动的频率并于输入的正弦曲线频率比较。

（4）根据实验情况，改变输入的正弦曲线频率的大小，重复上述，纪录实验数据。

（5）在仿真代码基础上，实现正弦激励代码，计算出实验结果对应的理论结果。对比分析理论和实验结果的差异。完成实验报告。

实验二：双自由度系统数学建模及仿真 1 实验目的

（1）熟悉双自由度质量-弹簧-阻尼系统并进行数学建模； （2）了解MATLAB 软件编程，学习编写系统的仿真代码； （3）进行双自由度系统的仿真动态响应分析。 2 实验原理 2.1 数学建模

双自由度质量-弹簧-阻尼系统，如上图所示。由两个质量为m 1和m 2的滑块、

两个刚度系数为k 1和k 2的弹簧和两个阻尼系数为c 1和c 2的阻尼器组成。系统输入：作用在滑块上的力f (t )。系统输出：滑块的位移x 1(t )和x 2(t )。 建立力学平衡方程：

112111111222122211122211()

m x c x c x k x k x f t m x c x c x c x k x k x k x ??

??

?????

?+-+-=???++-++-=? 2.2 固有频率

将动力学方程写成矩阵形式：

1111

1111211211222201()0()()0x m c c k k x x f t m c c c k k k x x x ??????????--???????????????+?+?=???????????????-+-+???

???????????

得到系统的质量矩阵M 和刚度矩阵K 。 解行列式可得固有频率方程：2||0K M ω-= 可计算出固有频率方程：

2

1121,2

120.5[k k k m m ω+=+

两个振动模态，两个固有频率：高模和低模。 2.3 解耦

通过数学变换将微分方程变化为以下形式：

2

1111111222222222()

2()

y y y f t y y y f t ζωω?ζωω????

???

?++=???++=? 注意：y 1和y 2不是滑块的位移。滑块的位移x 1(t )和x 2(t )是y 1和y 2的函数。

3 动力学仿真

根据数学模型，使用龙格库塔方法ODE45求解，任意输入下响应结果。 仿真代码见附件 4 实验

4.1 固有频率分析

（1）将实验台设置为双自由度质量-弹簧-阻尼系统，第一个滑块没有阻尼器可以不接，认为第一个阻尼为零。

（2）闭合控制器开关，点击setup 菜单，选择Control Algorithm ，设置选择Continuous Time Control ，Ts=0.0042,然后OK 。点击Command 菜单，选择Trajectory ，选取step,进入set-up,选取Open Loop Step 设置(0)counts, dwell time=3000ms, (1)rep, 然后OK 。此步是为了使控制器得到一段时间的数据，并不会驱动电机运动。

（2）点击Data 菜单，选择Data Acquisition,设置分别选取Encoder#1，Encoder#2，然后OK 离开；从Utility 菜单中选择Zero Position 使编码器归零。 （4）从Command 菜单中选择Execute ，用手将质量块1移动到2.5cm 左右的位置（注意不要使质量块碰触移动限位开关），点击Run, 大约1秒后，放开手使其自由振荡，在数据上传后点击OK 。

（5）点击Plotting 菜单，选择Setup Plot ，分别选取Encoder #1 Position ，Encoder #2

Position ；然后点击Plotting 菜单，选择Plot Data ，则将显示质量块1，2的自由振动响应曲线。 （6）实验数据纪录：

实验条件：滑块质量m 1和m 2，弹簧刚度k 1和k 2，阻尼系数c 1和c 2。

实验数据：时间-滑块1位移数据；时间-滑块2位移数据。

问题1：两个滑块位移的频率测量值是高模和低模频率么？

问题2：实际的机械系统是多自由度的，如何通过实验法测试系统固有频率？（7）实验报告。关键点是理论和实验结果对比分析。

4.2幅频特性实验

（1）点击Command菜单，选择Trajectory，选取Sinuscidal,进入set-up,选取Open Loop Step设置(200counts)Amplitude, Frequency(2Hz),Repetition（8），然后OK。

（2）从Utility菜单中选择Zero Position使编码器归零。从Command菜单中选择Execute，点击Run,在数据上传后点击OK。

（3）然后点击Plotting菜单，选择Plot Data，则将显示质量块1，2的受迫振动响应曲线。在响应曲线图上，即可测量出振动振幅。

问题1：单自由度和双自由度系统的幅频特性有何差异？

问题2：高模贡献分析。

实验三：PID控制

1 实验目的

（1）学习PID闭环控制结构和系统闭环传递函数计算；

（2）PID控制器参数设计；

（3）控制性能分析。

2 实验原理

上图给出闭环控制系统原理框图。单自由度质量-弹簧-阻尼系统结构下，断开弹簧和阻尼，仅仅保留滑块质量m。

m x f ??

=

电控箱可以看做比例增益k hw 。

hw f k u =

其中：u 是控制器输出。PID 控制：

p i d k e k edt k e u ?

++=?

其中：e 是比较器输出，参考输入与实际输出的偏差值。

e r x =-

根据全部上式，可得闭环结构微分方程：

hw d hw p hw i hw d hw p hw i m x k k x k k x k k x k k r k k r k k r ???

??

?

??

?

+++=++

对应的传递函数：

232()()()()

()()()()

hw d hw p hw i hw d hw p hw i k k s k k s k k x s r s ms k k s k k s k k ++=

+++ 3 PID 设计

PID 控制器中设置积分因子k i 为零，则为PD 控制。传递函数变为：

2()()

hw p hw d

hw p

hw d k k k k s x s m m k k k k r s s s m m

+

=++

闭环特征方程是分母：

ω=

2hw d k k m ζω=

=

设计频率

）的控制器设计。

4 实验

4.1 频率

（1）在控制器断开的情况下，拆除与质量块1连接的弹簧，使其余元件远离质量块1的运动范围，为其安装4个500g的铜块，加上小车本身的质量，标定总质量

m=2.6kg。

（2）实验标定k hw值：根据估计出的k hw值，设置控制器k i=0和k d=0，调整k p来估计系统系统频率ω=4Hz。注意： k p不能大于0.08。

（3）闭合控制器开关；点击Data菜单，选择Data Acquisition,设置选取Encoder#1和Commanded Position information；点击Command菜单，选择Trajectory，选取step，设置(0)counts,dwell time=3000ms,(1)rep。

（4）点击setup菜单，选择Control Algorithm，设置选择Continuous Time Control，Ts=0.0042,选取PID，进入Setup Algorithm，输入k p的值（k i=0和k d=0）（输入的值不能大于k p=0.08），然后OK。移动质量块1到-0.5cm的位置（规定，朝电机方向为负）选择Implement Algorithm，然后Ok。

注意：从此步开始的每一步，要进行下一步之前，都要与运动装置保持一定的安全距离；选择Implement Algorithm后，控制器将会立即加载，若出现不稳定的或者很大的控制信号时，运动装置可能反应很剧烈；若加载后，系统看上去稳定，要先用一轻质不尖锐的物体轻轻碰触质量块以验证其稳定性。

（5）点击Command菜单，进入Execute，用手将质量块移动到2cm左右的位置，点击Run,移动质量块大约到3cm位置，然后释放（不要拿着质量块多于一秒，以免电机过热而断开控制）。

（6）点击Plotting菜单，选择Setup Plot，选取Encoder #1；然后点击Plotting 菜单，选择Plot Data，则将显示质量块1的时间响应曲线。

试想一下，若将比例增益系数k p增加一倍，则系统的响应频率将有什么变化？

4.2 阻尼

（1）确定k d的值（不能大于0.04），使得k hw k d=50N/(m/s)，重复步骤4，除了输入kd的值和k i=0和k p=0。

（2）先用尺子检查系统的稳定性，然后用手来回的移动质量块来感受系数k d带来的粘性阻尼的影响（注意不要极度的迫使质量块运动，以免电机过热而断开控制）。

（3）增大的k d值（k d<0.04），重复以上步骤，看能否感受到阻尼的增大？4.3 位置控制

（1）控制电机驱动滑块1移动到某一个特定的位置，获取时间-位移数据后，计算闭环控制特性，包括：超调量、上升时间、调节时间、稳态误差。

（2）使用仿真方法计算得到同样条件下的动态响应结果，对比分析理论与实验结果的差异。

4.4 正弦激励

（1）使用正弦输入驱动滑块作正弦运动，获取时间-位移数据后，计算闭环控制系统幅频特性。

（2）使用仿真方法计算得到同样条件下的幅频响应结果，对比分析理论与实验结果的差异。

4.5 柔性结构控制

（1）实验台只有质量块，为刚性结构机械系统。如果安装上一根弹簧，则组成简单的柔性结构系统（具有柔性环节）。此时机械系统由于能量交互可能产生振动。如何设计控制器来动态控制机械振动。

实验四：扰动抑制

1 实验目的

（1）外扰动设置；

（2）PID控制器设计；

（3）控制性能分析。

2 实验原理

上图给出闭环控制系统原理框图。单自由度质量-弹簧-阻尼系统结构下，断开弹簧和阻尼，仅仅保留滑块质量m。位移输出受到外扰动d干扰。

m x f d ??

=+

电控箱可以看做比例增益k hw 。

hw f k u =

其中：u 是控制器输出。PID 控制：

p i d k e k edt k e u ?

++=?

其中：e 是比较器输出，参考输入与实际输出的偏差值。

e r x =-

根据全部上式，可得闭环结构微分方程：

hw d hw p hw i hw d hw p hw i m x k k x k k x k k x k k r k k r k k r d ???

??

?

??

?

?

+++=+++

对应的传递函数：

23232()()()()()()

()()()

()()()

hw d hw p hw i hw d hw p hw i hw d hw p hw i k k s k k s k k s

x s r s d s ms k k s k k s k k ms k k s k k s k k ++=+

++++++

3 PID 设计

PID 控制器中设置积分因子k i 为零，则为PD 控制。传递函数变为：

2()()

hw p

hw d

hw p

hw d k k k k s x s m m k k k k r s s s m m

+=++

闭环特征方程是分母：

ω=

2hw d k k m ζω=

=

附件1: 实验1仿真代码

% exp 1

clc

clear

global SK SC SM SF;

sampling=1/1000;

SK=200;

SC=0;

SM=2.6;

State=zeros(1,2);

State(1)=1;

State(2)=0.0;

sss=zeros(1,1);

for k=1:fix(5/sampling)

t=k*sampling;

sss(k,1)=t; %time

SF=0;

sss(k,2)=SF;

[TimeOdeArray,VarOdeArray] = ode45( @mdlDerivatives1, [t t+sampling],State); [m,n]=size(TimeOdeArray);

TimeAtEndOfArray = TimeOdeArray(m,1);

if( abs(TimeAtEndOfArray - (t+sampling) ) >= abs(0.001*sampling) )

warning('numerical integration failed'); break; end

VAR = VarOdeArray(m, : );%ò???2?3¤?áê?μ?×′ì?

State=VAR;%3?ê?×′ì??üD?￡???DD??ò?′?μü′ú

sss(k,3)=VAR(1); % response

sss(k,4)=VAR(2); % velocity

end

plot(sss(:,1),sss(:,2),'r',sss(:,1),sss(:,3),'k',sss(:,1),sss(:,4),'b');

% mdlDerivatives1

function dx=mdlDerivatives1(T,x)

global SK SC SM SF;

% y(1)=x;

% y(2)=xd;

dx=zeros(2,1);

%sloshing dynamics

dx(1)=x(2);

dx(2)=SF/SM-SK/SM*x(1)-SC/SM*x(2);

附件2: 实验2仿真代码

% exp 2

clc

clear

global SK1 SK2 SC1 SC2 SM1 SM2 SF;

sampling=1/1000;

SK1=200;

SK2=200;

SC1=0;

SC2=0;

SM1=2.6;

SM2=2.6;

State=zeros(1,4);

State(1)=1.0;

State(3)=0.0;

sss=zeros(1,1);

for k=1:fix(5/sampling)

t=k*sampling;

sss(k,1)=t; %time

SF=0;

sss(k,2)=SF;

[TimeOdeArray,VarOdeArray] = ode45( @mdlDerivatives2, [t t+sampling],State); [m,n]=size(TimeOdeArray);

TimeAtEndOfArray = TimeOdeArray(m,1);

if( abs(TimeAtEndOfArray - (t+sampling) ) >= abs(0.001*sampling) )

warning('numerical integration failed'); break; end

VAR = VarOdeArray(m, : );%ò???2?3¤?áê?μ?×′ì?

State=VAR;%3?ê?×′ì??üD?￡???DD??ò?′?μü′ú

sss(k,3)=VAR(1); % response 1

sss(k,4)=VAR(2); % velocity 1

sss(k,5)=VAR(3); % response 2

sss(k,6)=VAR(4); % velocity 2

end

plot(sss(:,1),sss(:,2),'r',sss(:,1),sss(:,3),'k',sss(:,1),sss(:,5)-sss(:,3),'b');

% mdlDerivatives2

function dx=mdlDerivatives2(T,x)

global SK1 SK2 SC1 SC2 SM1 SM2 SF;

% x(1)=x1;

% x(2)=x1d;

% x(3)=x2;

% x(4)=x2d;

dx=zeros(4,1);

%sloshing dynamics

dx(1)=x(2);

dx(2)=(SF-SC1*x(2)+SC1*x(4)-SK1*x(1)+SK1*x(3))/SM1;

dx(3)=x(4);

dx(4)=(-1*SC2*x(4)-SC1*x(4)+SC1*x(2)-SK1*x(3)-SK2*x(3)+SK1*x(1))/SM2;

附件3 快速傅里叶变换FFT计算频率代码%frequency response

sampling=1/1000;

Fs=fix(1/sampling);

data=sss(1:fix(5/sampling),3);

m=length(data);

nfft=2^nextpow2(m);

y=fft(data,nfft);%?μà?ò?±???Ayy=abs(y)*2/nfft;

f0=(0:nfft/2-1)*Fs/nfft; % Hz , 1/s

f1=f0*2*pi; % rad/s

zzz(:,1)=f0;

zzz(:,2)=f1;

zzz(:,3)=Ayy(1:nfft/2);

plot(zzz(:,1),zzz(:,3));

xlabel('frequency/Hz');

ylabel('magnitude');

title('freqency response');

1 Command/Trajectory

输入命令

2 Command/Execute

执行

3 plotting/Setupplot

绘图运动轨迹

4 Setup改变控制器