汽车振动学实验报告

《汽车振动学》实验报告姓名: 舒鹏

学号: 20110802

班级: 2011级交通运输2班

指导教师: 蒋淑霞

实验主题:基于Matlab的汽车振动分析

实验一：单自由度阻尼振动分析：包括欠阻尼 、临界阻尼

和过阻尼

一、 实验内容

（1）将老师给的源程序读懂（包含临界阻尼和欠阻尼）。

（2）参考上述程序，改变阻尼比

和画出不同阻尼比下的响应曲线共10条，并比较各响应曲线特性。 （3）同时考虑过阻尼曲线。

二、 实验前期准备与内容分析 阻尼自由振动：

前面讲的自由振动系统没有考虑阻力的影响，实际由于阻力的存在系统机械能不可能守恒，振动中这种阻力称为阻尼，如摩擦阻尼、电磁阻尼、介质阻尼和结构阻尼。对于实际系统中的阻尼很难确定，工程中最常见的一种阻尼力学模型是粘性阻尼，如在流体中运动、在润滑表面的滑动。 粘性阻尼力与相对速度成正比:

c ：粘性阻尼系数，或阻尼系数，单位：N ·s/m

建立平衡位置坐标系，受力分析：阻尼力、弹性力、重力.

得到动力学方程:

其中，固有频率： 相对阻尼系数:

三、 实验代码与结果 第一种情况：零界阻尼

ksai=1 %零界阻尼,ksai 表示为阻尼比 k=10^3; m=4;

w0=sqrt(k/m); %固有频率 v0=10; %初始速度 x0=5; %初始位移 t=(0:300)/100;

x1=exp(-w0*t).*(x0+v0*t) axes(handles.axes1)

plot(t,x1) %在图形axes1中输出图像 xlabel('时间 （s ）'); %x 坐标表示时间 ylabel('位移 （mm ）'); %y 坐标表示位移

第二种情况：欠阻尼

ksai=0.2 %阻尼比为0.2

k=10^3;

1,...,2.0,1.0=

ε

m=4;

x0=5 %初始位移

v0=10 %初始速度

w0=sqrt(k/m); %固有频率

wd=sqrt(1-ksai^2)*w0; %阻尼固有频率

t=(0:300)/100;

x2=exp(-ksai*w0*t).*(x0*cos(wd*t)+((v0+ksai*w0*x0)/wd)*sin(wd*t ))

y=x0*exp(-ksai*w0*t) %渐近线

axes(handles.axes2)

plot(t,x2,'k',t,y,'g--') ; 在axes2输出图像

hold on;

ksai=0.4 %阻尼比0.4

k=10^3;

m=4;

x0=5 %初始位移

v0=10 %初始速度

w0=sqrt(k/m); %固有频率

wd=sqrt(1-ksai^2)*w0; %阻尼固有频率

t=(0:300)/100;

x2=exp(-ksai*w0*t).*(x0*cos(wd*t)+v0*sin(wd*t))

axes(handles.axes2)

plot(t,x2,'y')

hold on;

ksai=0.6 %阻尼比0.6

k=10^3;

m=4;

x0=5 %初始位移

v0=10 %初始速度

w0=sqrt(k/m); %固有频率

wd=sqrt(1-ksai^2)*w0; %阻尼固有频率

t=(0:300)/100;

x2=exp(-ksai*w0*t).*(x0*cos(wd*t)+v0*sin(wd*t))

axes(handles.axes2)

plot(t,x2,'m')

hold on;

legend('欠阻尼0.2','零界线','欠阻尼0.4','欠阻尼0.6') %对曲线进行命名

xlabel('时间（s）');

ylabel('位移（mm）');

第三种情况：过阻尼

ksai=2 %阻尼比为2

k=10^3;

m=4;

x0=5 %初始位移

v0=10 %初始速度

w0=sqrt(k/m); %固有频率

w1=w0*sqrt(ksai^2-1);

t=(0:300)/100

shw1t=(exp(w1*t)-exp(-w1*t))/2

chw1t=(exp(w1*t)+exp(-w1*t))/2

x3=exp(-ksai*w0*t).*(x0*chw1t+((v0+ksai*w0*x0)/w1)*shw1t) axes(handles.axes3);

plot(t,x3)

xlabel('时间（s）');

ylabel('位移（mm）');

在gui中输出结果：

实验二：单自由度简谐振动幅频相频特性曲线

一、实验内容

1.观察单自由度简谐振动的幅频相频特性曲线。

2.掌握放大因子β随着频率比变化的规律及曲线随着阻尼比ζ变化的规律。

二、实验前期准备与内容分析

稳态响应特性

以s为横坐标，画出幅频特性曲线：

三、实验代码与结果

幅频曲线：

s=(0:300)/100; %频率比

ksai=0.1 %阻尼比0.1

b1=1./sqrt((1-s.^2).^2+(2*ksai*s).^2)

plot(s,b1); % 输出图像

hold on;

ksai=02 %阻尼比0.2

b1=1./sqrt((1-s.^2).^2+(2*ksai*s).^2)

plot(s,b1,'k--')

hold on;

ksai=0.3 %阻尼比0.3

b1=1./sqrt((1-s.^2).^2+(2*ksai*s).^2)

plot(s,b1,'b')

hold on;

ksai=0.4 %阻尼比0.4

b1=1./sqrt((1-s.^2).^2+(2*ksai*s).^2)

plot(s,b1,'y')

hold on;

ksai=0.5 %阻尼比0.5

b1=1./sqrt((1-s.^2).^2+(2*ksai*s).^2)

plot(s,b1,'m')

hold on;

b1=(0:60)/10

plot(1,b1)

legend('ksai=0.1','ksai=0.2','ksai=0.3','ksai=0.4','ksai=0.5') xlabel('间歇力频率与固有频率之比');

ylabel('响应振幅与静位移之比 ');

相频曲线：

s=(0:500)/100 ; %频率比

ksai=0.1 ; %阻尼比0.3

A=atan(2*ksai*s./(1-s.^2));

for i=1:size(A,2);

if A(1,i)<0;

A(1,i)=A(1,i)+pi;

end

end

plot(s,A*180/pi,'k'); hold on;

ksai=0.2 ; %阻尼比0.2

A=atan(2*ksai*s./(1-s.^2));

for i=1:size(A,2);

if A(1,i)<0;

A(1,i)=A(1,i)+pi;

end

end

plot(s,A*180/pi,'m'); hold on;

ksai=0.3 ; %阻尼比0.3

A=atan(2*ksai*s./(1-s.^2));

for i=1:size(A,2);

if A(1,i)<0;

A(1,i)=A(1,i)+pi;

end

end

plot(s,A*180/pi,'b'); hold on;

ksai=0.4 ; %阻尼比0.4

A=atan(2*ksai*s./(1-s.^2));

for i=1:size(A,2);

if A(1,i)<0;

A(1,i)=A(1,i)+pi;

end

end

plot(s,A*180/pi,'y'); hold on;

ksai=0.5 %阻尼比0.5

A=atan(2*ksai*s./(1-s.^2));

for i=1:size(A,2);

if A(1,i)<0;

A(1,i)=A(1,i)+pi;

end

end

plot(s,A*180/pi); hold on;

legend('阻尼比0.1','阻尼比0.2','阻尼比0.3','阻尼比0.4','阻尼比0.5')

xlabel('频率比');

ylabel('相位角');

幅频曲线图

相频曲线图：

实验三 基于MatLab 的偏心质量引起的强迫振动分析

一、实验内容

? 基于MatLab 的偏心质量引起的强迫振动分析。 ? （1）参考书上的内容p6-p8，特别是公式 ? 1.2.16和1.2.17，求动态放大因子 随频率 ? 比

n ωω 和阻尼比 km

2c =ε 的变化曲线。 ? （2）求响应相位角滞后于激振力的相位角 ? 随频率比

n ωω 和阻尼比 km

2c

=ε 的变化曲线。 二、实验前期准备与内容分析

旋转机械如电动机、水泵、汽车发动机等，由于偏心质量而引起强迫振动是很普遍的现象。

偏心质量引起的强迫振动分析

()()

()2222

022222212

214ζλ

λλωωω+-=+-

=B n p m e m B

2

2

2ω

ω

ψ-=p n arctg

基于MatLab的偏心质量引起的强迫振动分析

四、实验代码与结果

实验二程序与实验一程序相同，只是只是代表物理含义名称不同s=(0:300)/100; %频率比

ksai=0.1 %阻尼比0.1

b1=1./sqrt((1-s.^2).^2+(2*ksai*s).^2)

plot(s,b1); % 输出图像

hold on;

ksai=02 %阻尼比0.2

b1=1./sqrt((1-s.^2).^2+(2*ksai*s).^2)

plot(s,b1,'k--')

hold on;

ksai=0.3 %阻尼比0.3

b1=1./sqrt((1-s.^2).^2+(2*ksai*s).^2)

plot(s,b1,'b')

hold on;

ksai=0.4 %阻尼比0.4

b1=1./sqrt((1-s.^2).^2+(2*ksai*s).^2)

plot(s,b1,'y')

hold on;

ksai=0.5 %阻尼比0.5

b1=1./sqrt((1-s.^2).^2+(2*ksai*s).^2)

plot(s,b1,'m')

hold on;

b1=(0:60)/10

plot(1,b1)

legend('ksai=0.1','ksai=0.2','ksai=0.3','ksai=0.4','ksai=0.5') xlabel('频率比');

ylabel('放大因子 ');

s=(0:500)/100 ; %频率比

ksai=0.1 ; %阻尼比0.3

A=atan(2*ksai*s./(1-s.^2));

for i=1:size(A,2);

if A(1,i)<0;

A(1,i)=A(1,i)+pi;

end

end

plot(s,A*180/pi,'k'); hold on;

ksai=0.2 ; %阻尼比0.2

A=atan(2*ksai*s./(1-s.^2));

for i=1:size(A,2);

if A(1,i)<0;

A(1,i)=A(1,i)+pi;

end

end

plot(s,A*180/pi,'m'); hold on;

ksai=0.3 ; %阻尼比0.3

A=atan(2*ksai*s./(1-s.^2));

for i=1:size(A,2);

if A(1,i)<0;

A(1,i)=A(1,i)+pi;

end

end

plot(s,A*180/pi,'b'); hold on;

ksai=0.4 ; %阻尼比0.4

A=atan(2*ksai*s./(1-s.^2));

for i=1:size(A,2);

if A(1,i)<0;

A(1,i)=A(1,i)+pi;

end

end

plot(s,A*180/pi,'y'); hold on;

ksai=0.5 %阻尼比0.5

A=atan(2*ksai*s./(1-s.^2));

for i=1:size(A,2);

if A(1,i)<0;

A(1,i)=A(1,i)+pi;

end

end

plot(s,A*180/pi); hold on;

legend('阻尼比0.1','阻尼比0.2','阻尼比0.3','阻尼比0.4','阻尼比

0.5')

xlabel('频率比');

ylabel('相位角');

幅频曲线图

实验四：基于Matlab 的三自由度振动分析

一、实验内容

?参考源程序1 motaiqiujie.m

?符号方程的求解，解特征方程得到三个振动固有频率；

?源程序2 chapter2_3freed.m

?（1）首先读懂程序

?（2）系统M和K矩阵如下，且m=1000kg；k=10^6，求固有频率和三个振型的实频和虚频曲线、幅频和相频曲线和Nyquist图。理解主振动含义。

二、实验前期准备与内容分析

Matlab求特征值和特征向量:

三、实验代码与结果

a=0.001;b=0.001;

M=[100,0,0;0,100,0;0,0,100];

K1=[2,-1,0;-1,2,-1;0,-1,2];

K=1.0e+5*K1;

C=a*M+b*K; %比例阻尼

[V,D]=eig(K,M);

% [V,D] = EIG(A,B) produces a diagonal matrix D of generalized

% eigenvalues and a full matrix V whose columns are the

% corresponding eigenvectors so that A*V = B*V*D.

for j=1:1:3

omega(j)=sqrt(D(j,j));

for k=1:1:3

Y(k,j)=V(k,j)/V(3,j);

end

end

omega

Y

w=0:0.1:70; %频率取值范围和步长H=zeros(701,3,3);

for j=1:1:701

Z=K-w(j)^2*M+w(j)*C*i; %动刚度矩阵

HH=inv(Z); %频响函数矩阵

H(j,:,:)=HH;

end;

Hr=real(H);Hi=imag(H);

Habs=abs(H);Hang=(180/pi)*angle(H);

%下面分别是三个振型的实频和虚频曲线；

figure(1)

subplot(3,2,1),plot(w,Hr(:,1,1)),grid on

title('实频曲线1'),xlabel('\omega'),ylabel('Hr')

subplot(3,2,2),plot(w,Hi(:,1,1)),grid on

title('虚频曲线1'),xlabel('\omega'),ylabel('Hi')

subplot(3,2,3),plot(w,Habs(:,1,1)),grid on

title('幅频曲线1'),xlabel('\omega'),ylabel('|H|')

subplot(3,2,4),plot(w,Hang(:,1,1)),grid on

title('相频曲线1'),xlabel('\omega'),ylabel('\psi')

subplot(3,2,6);

plot(Hr(:,1,1),Hi(:,1,1)),grid on

title('乃奎斯特圆1'),xlabel('Re(H)'),ylabel('Im(H)')

%画跨点的频响

figure(2)

subplot(3,2,1);plot(w,Hr(:,2,1)),grid on

title('实频曲线2'),xlabel('\omega'),ylabel('Hr')

subplot(3,2,2);plot(w,Hi(:,2,1)),grid on

title('虚频曲线2'),xlabel('\omega'),ylabel('Hi')

subplot(3,2,3);plot(w,Habs(:,2,1)),grid on

title('幅频曲线2'),xlabel('\omega'),ylabel('|H|')

subplot(3,2,4);plot(w,Hang(:,2,1)),grid on

title('相频曲线2'),xlabel('\omega'),ylabel('\psi')

subplot(3,2,6);

plot(Hr(:,2,1),Hi(:,2,1)),grid on

title('乃奎斯特圆2'),xlabel('Re(H)'),ylabel('Im(H)')

%画跨点的频响

figure(3)

subplot(3,2,1);plot(w,Hr(:,3,1)),grid on

title('实频曲线3'),xlabel('\omega'),ylabel('Hr')

subplot(3,2,2);plot(w,Hi(:,3,1)),grid on

title('虚频曲线3'),xlabel('\omega'),ylabel('Hi')

subplot(3,2,3);plot(w,Habs(:,3,1)),grid on

title('幅频曲线3'),xlabel('\omega'),ylabel('|H|')

subplot(3,2,4);plot(w,Hang(:,3,1)),grid on

title('相频曲线3'),xlabel('\omega'),ylabel('\psi')

subplot(3,2,6);

plot(Hr(:,3,1),Hi(:,3,1)),grid on

title('乃奎斯特圆3'),xlabel('Re(H)'),ylabel('Im(H)')

w=sqrt(eig(K/M));

实验五：多自由度汽车振

动分析

一、实验内容

给定发动机的参数书上的题3.5（p98）未求和答案，请优化程序并计算结果。

二、实验前期准备与内容分析

某4支撑发动机悬置系统，已知如下参数：

质量m 转动惯量J 刚度矩阵K 线位移矩阵X 角位移矩阵θ

求系统的质量矩阵、刚度矩阵、固有频率和振型，并指出系统的耦合情况。

三、实验代码与结果

clear all; m=0.51;

Jx=236;Jy=1087;Jz=957;Jxy=0;Jyz=0;Jxz=0;

for i=1:4

ou(i)=0; ov(i)=0; os(i)=0;

end

ku(1)=249; kv(1)=249; ks(1)=878;

ku(2)=249; kv(2)=249; ks(2)=878;

ku(3)=1849; kv(3)=623; ks(3)=971;

ku(4)=1849; kv(4)=623; ks(4)=971;

x(1)=53.17*0.01; y(1)=26.18*0.01; z(1)=-14.6*0.01;

x(2)=53.17*0.01; y(2)=-27.42*0.01; z(2)=-14.6*0.01;

x(3)=-27.33*0.01; y(3)=29.38*0.01; z(3)=-11.6*0.01;

x(4)=-27.33*0.01; y(4)=-30.62*0.01; z(4)=-11.6*0.01;

ov(3)=-3*pi/180; ov(4)=-3*pi/180; ou(1)=-40*pi/180;

ou(2)=40*pi/180; ou(3)=5*pi/180; ou(4)=5*pi/180;

%求系统的质量矩阵、刚度矩阵、固有频率和振型

%求解过程与答案

for i=1:6

for j=1:6 M(i,j)=0;

end

end

M(1,1)=m; M(2,2)=m; M(3,3)=m;

M(4,4)=Jx; M(5,5)=Jy; M(6,6)=Jz;

M(4,5)=-Jxy; M(5,4)=-Jxy; M(4,6)=-Jxz;

M(6,4)=-Jyz; M(5,6)=-Jyz; M(6,5)=-Jyz;M

sum1=zeros(6,6) ;

for k=1:4

Cu(:,:,k)=[1 0 0;0 cos(ou(k)) sin(ou(k));0 -sin(ou(k)) cos(ou(k))];

Cv(:,:,k)=[cos(ov(k)) 0 -sin(ov(k)); 0 1 0; sin(ov(k)) 0 cos(ov(k))]; Cs(:,:,k)=[cos(os(k)) sin(os(k)) 0 ; -sin(os(k)) cos(os(k)) 0;0 0 1]; C(:,:,k)=Cs(:,:,k)*Cu(:,:,k)*Cv(:,:,k);

D(:,:,k)=[ku(k) 0 0;0 kv(k) 0; 0 0 ks(k)]

T(:,:,k)=[1 0 0 0 z(k) -y(k);0 1 0 -z(k) 0 x(k); 0 0 1 y(k) -x(k) 0];

K(:,:,k)=(T(:,:,k))'*(C(:,:,k))'*D(:,:,k)*C(:,:,k)*T(:,:,k); end

KK=sum(K,3)

C=inv(M)*KK

[A,B]=eig(C) %A表示主振型；B表示固

有频率平方组成对角阵

w=diag(sqrt(B)) %固有频率C=K*inv(M) f1=w/2/pi

其中，质量矩阵为M、刚度矩阵为KK、固有频率W.

实验六：利用行列式值求矩阵特征根！

一、实验内容

利用行列式值求矩阵特征根！

二、实验前期准备与内容分析

无阻尼三自由度系统运动微分方程：

m、k、x分别为气门、挺杆和摇臂等效质量、

等效刚度和垂直位移。

F0-气门推力

--摇臂转角

系统运动微分方程求解：

特征频率方程：

给出：

求出固有频率：

三、实验代码与结果

%z,x,y方向的传感器标定值分别为10.19，10.22，10.74 Fs=1000;

load y

data(:,2)=y(2001:4000);

data(:,3)=y(4001:6000);

data(:,4)=y(6001:8000);

x=data(:,2)/10.22;

y=data(:,3)/10.74;

z=data(:,4)/10.19;

figure(1)

subplot(3,1,1)

plot(x)

title('x')

subplot(3,1,2)

plot(y)

title('y')

subplot(3,1,3)

plot(z)

title('z')

figure(2)

subplot(3,1,1);

[pz,fz] = filtbank(z,Fs);

hz = bankdisp(pz,fz);

title('座椅Z方向1/3倍频程分析');

wz=[0.63 0.71 0.8 0.9 1 1 1 1 0.8 0.63 0.5 0.4 0.315 0.25 0.2 0.16 0.125 0.1];

subplot(3,1,2);

[px,fx] = filtbank(x,Fs);

hx = bankdisp(px,fx);

title('座椅X方向1/3倍频程分析');

wx=[1 1 0.8 0.63 0.50 0.40 0.315 0.25 0.20 0.16 0.125 0.10 0.08 0.063 0.05 0.04 0.0315 0.025];

subplot(3,1,3);

[py,fy] = filtbank(y,Fs);

hy = bankdisp(py,fy);

title('座椅Y方向1/3倍频程分析');

wy=[1 1 0.8 0.63 0.50 0.40 0.315 0.25 0.20 0.16 0.125 0.10 0.08 0.063 0.05 0.04 0.0315 0.025];

wz2=wz.^2;

pz2=pz.^2;

az2=wz2.*pz2;

az2sum=sum(az2);

az=sqrt(az2sum);

wx2=wx.^2;

px2=px.^2;

ax2=wx2.*px2;

ax2sum=sum(ax2);

ax=sqrt(ax2sum);

wy2=wy.^2;

py2=py.^2;

ay2=wy2.*py2;

ay2sum=sum(ay2);

ay=sqrt(ay2sum);

ao=10^-6;

awo=((1.4*ax)^2+(1.4*ay)^2+az^2)^(1/2)

leq=20*log10(awo/ao)

leqx=20*log10(ax/ao)

leqy=20*log10(ay/ao)

leqz=20*log10(ay/ao)

运行计算结果

leq = 62.3969 leqx = 53.2024

leqy = 50.7895 leqz =50.7895

awo = 0.0013 ax = 4.5721e-004

ay =3.4632e-004 az=0.0010

实验七扭转振动分析的传递矩阵法一、实验内容

求三圆盘扭振系统的固有频率和扭振模态

二、实验前期准备与内容分析

三、实验代码与结果

clear all;

syms wn;

%A(2,2)=1-(wn).^2/10^4;A

A=[1 1/(2*10^7);(-2*10^3)*wn^2 1-wn^2/10^4] B=[1 1*10^-7;-(1)*10^3*(wn)^2 1-(wn)^2/10^4] C=[1; -500*wn^2]

M=A*B*C