振动理论及应用
大作业
题目：基于MATLAB的三自由度
弹簧质量系统的振动分析
姓名：
专业：
学号：
一、题目：
已知图示的三自由度弹簧质量系统，试编写MATLAB程序：
（1）求固有频率、主振型及正则振型；
（2）对初始条件的自由响应；
（3）对外激励的稳态响应。
参数选取：（1）各个质量值（2）各段刚度值（3）初始条件（4）简谐激励。
二、求解过程和MATLAB程序
取k1=k2=k3=k4=k，m1
=
m2=m，m3=2m；初始条件为，简谐激励为。
1、求解固有频率、主振型及正则振型
建立名为work.m的m文件，并输入以下命令：
%定义刚度矩阵和质量矩阵
k=1;m=1;
k1=k;k2=k;k3=k;k4=k;
m1=m;m2=m;m3=2*m;
k11=k1+k2;k12=-k2;k13=0;
k21=-k2;k22=k2+k3;k23=-k3;
k31=0;k32=-k3;k33=k3;
K=[k11,k12,k13;k21,k22,k23;k31,k32,k33];
M=[m1,0,0;0,m2,0;0,0,m3];
%求特征值和特征向量
R=inv(M)*K;
D1=eig(R);
[p2,d]=sort(D1);
[V,Dm]=eig(R);
%固有频率
for
i=1:3
p(i,1)=sqrt(p2(i,1));
end
%求主振型和正则振型
for
i=1:3
for
j=1:3
A1(i,j)=V(i,d(j))/V(1,d(j));
end
end
Ap=A1;
Mp=Ap＇*M*Ap;
Kp=Ap＇*K*Ap;
for
i=1:3
AN(:,i)=Ap(:,i)/sqrt(Mp(i,i));
end
MN=AN＇*M*AN;
KN=AN＇*K*AN;
p
Ap
AN
运行得到固有频率p，主振型矩阵Ap，正则振型矩阵AN如下：
p
=
0.3560
1.1281
1.7609
Ap
=
1.0000
1.0000
1.0000
1.8733
0.7275
-1.1007
2.5092
-0.4708
0.2116
AN
=
0.2418
0.7120
0.6592
0.4530
0.5180
-0.7256
0.6068
-0.3352
0.13952、求解对初始条件的自由响应
（1）建立函数文件FunFree.m，并输入如下命令：
function
dy=FunFree(t,y)
k=1;m=1;
k1=k;k2=k;k3=k;k4=k;
m1=m;m2=m;m3=2*m;
k11=k1+k2;k12=-k2;k13=0;
k21=-k2;k22=k2+k3;k23=-k3;
k31=0;k32=-k3;k33=k3+k4;
dy=zeros(6,1);
dy(1)=y(4);
dy(2)=y(5);
dy(3)=y(6);
dy(4)=-1/m1*(k11*y(1)+k12*y(2)+k13*y(3));
dy(5)=-1/m2*(k21*y(1)+k22*y(2)+k23*y(3));
dy(6)=-1/m3*(k31*y(1)+k32*y(2)+k33*y(3));
（2）建立名为MainFree.m的m文件，并输入如下命令：
t=[0,50];
y0=[1;0;0;0;1;0];
[T,X]=ode45(@FunFree,t,y0);
figure;
subplot(3,1,1);
plot(T,X(:,1),＇r-＇);
subplot(3,1,2);
plot(T,X(:,2),＇b-＇);
subplot(3,1,3);
plot(T,X(:,3),＇k-＇);
运行得对初始条件的自由响应如下：
3、求解对外激励的稳态响应
建立函数文件Forced.m，并输入如下命令：
k=1;m=1;w=1;
k1=k;k2=k;k3=k;k4=k;
m1=m;m2=m;m3=2*m;
k11=k1+k2;k12=-k2;k13=0;
k21=-k2;k22=k2+k3;k23=-k3;
k31=0;k32=-k3;k33=k3+k4;
K=[k11,k12,k13;k21,k22,k23;k31,k32,k33];
M=[m1,0,0;0,m2,0;0,0,m3];
R=inv(M)*K;
D1=eig(R);
[p2,d]=sort(D1);
[V,Dm]=eig(R);
for
i=1:3
p(i,1)=sqrt(p2(i,1));
end
for
i=1:3
for
j=1:3
A1(i,j)=V(i,d(j))/V(1,d(j));
end
end
Ap=A1;
Mp=Ap＇*M*Ap;
Kp=Ap＇*K*Ap;
for
i=1:3
AN(:,i)=Ap(:,i)/sqrt(Mp(i,i));
end
MN=AN＇*M*AN;
KN=AN＇*K*AN;
F1=[1;0;0];
F2=[0;2;0];
F3=[0;0;3];
QN1=AN＇*F1*sin(w*t);
QN2=AN＇*F2*sin(2*w*t);
QN3=AN＇*F3*sin(3*w*t);
b1=[1/(p2(1,1)-w^2),0,0;
0,1/(p2(2,1)-w^2),0;
0,0,1/(p2(3,1)-w^2)];
b2=[1/(p2(1,1)-(2*w)^2),0,0;
0,1/(p2(2,1)-(2*w)^2),0;
0,0,1/(p2(3,1)-(2*w)^2)];
b3=[1/(p2(1,1)-(3*w)^2),0,0;
0,1/(p2(2,1)-(3*w)^2),0;
0,0,1/(p2(3,1)-(3*w)^2)];
XN1=b1*QN1;
XN2=b2*QN2;
XN3=b3*QN3;
T=20;
t=1:0.1:T;
X=AN*XN1+AN*XN2+AN*XN3;
figure;
subplot(3,1,1);
plot(t,X(1,:),＇r＇);
subplot(3,1,2);
plot(t,X(2,:),＇b＇);
subplot(3,1,3);
plot(t,X(3,:),＇k＇);
运行得对外激励的稳态响应如下：
点击咨询