动力学第三章.doc

第2章
function VTB2(m,c,k,x0,v0,tf,w,f0)
%单自由度系统得谐迫振动
clc
wn=sqrt(k/m);
z=c/2/m/wn;
lan=w/wn
wd=wn*sqrt(1-z^2);
A=sqrt(((v0+z*wn*x0)^2+(x0*wd)^2)/wd^2);
t=0:tf/1000:tf;
phi1=atan2(x0*wd,v0+z*wn*x0);
phi2=atan2(2*z*lan,1-lan^2);
B=wn^2*f0/k/sqrt((wn^2-w^2)^2+(2*z*wn*w)^2);
x=A*exp(-z*wn*t)、*sin(sqrt(1-z^2)*wn*t+phi1)+B*sin(w*t-phi2);
plot(t,x),grid
xlabel('时间(s)')
ylabel('位移')
title('位移与时间得关系')
function VTB1(m,c,k,x0,v0,tf)
%VTB1用来计算单自由度有阻尼自由振动系统得响应
%VTB1绘出单自由度有阻尼自由振动系统得响应图
%m为质量;c为阻尼;k为刚度;x0为初始位移;v0为初始速度;tf为仿真时间
%VTB1(zeta,w,x0,v0,tf)绘出单自由度有阻尼自由振动系统得响应图
%程序中z为阻尼系数ξ;wn为固有频率ωn;A为振动幅度;phi为初相位θ
clc
wn=sqrt(k/m);
z=c/2/m/wn;
wd=wn*sqrt(1-z^2);
fprintf('固有频率为%、3g、rad/s、\n',wd);
fprintf('阻尼系数为%、3g、\n',z);
fprintf('有阻尼得固有频率为%、3g、rad/s、\n',wd);
t=0:tf/1000:tf;
if z<1
A=sqrt(((v0+z*wn*x0)^2+(x0*wd)^2)/wd^2);
phi=atan2(x0*wd,v0+z*wn*x0)
x=A*exp(-z*wn*t)、*sin(wd*t+phi);
fprintf('A=%、3g\n',A);
elseif z==1
a1=x0;
a2=v0+wn*x0;
fprintf('a1=%、3g\n',a1);
fprintf('a2=%、3g\n',a2);
x=(a1+a2*t)、*exp(-wn*t);
else
a1=(-v0+(-z+sqrt(z^2-1))*wn*x0)/2/wn/sqrt(z^2-1);
a2=(v0+(z+sqrt(z^2-1))*wn*x0)/2/wn/sqrt(z^2-1);
fprintf('a1=%、3g\n',a1);
fprintf('a2=%、3g\n',a2);
x=exp(-z*wn*t)、*(a1*exp(-wn*sqrt(z^2-1)*t)+a2*exp(wn*sqrt(z^2-1)*t));
end
plot(t,x),grid
xlabel('时间(s)')
ylabel('位移')
title('位移与时间得关系')
function jzdd
%矩阵迭代法求系统得三阶固有频率与主阵型
clear all
close all
fid1=fopen('A-1','wt'); %建立主振型文件
fid2=fopen('B-1','wt'); %建立固有频率文件
%输入质量矩阵
M(1,1)=2;
M(2,2)=1、5;
M(3,3)=1;
%输入刚度矩阵
K(1,1)=5;K(1,2)=-2;K(2,1)=-2;K(2,2)=3;K(2,3)=-1;K(3,2)=-1;K(3,3)=1
%计算特征值与特征向量
D=inv(K)*M; %原始动力矩阵
A=ones(3,1); %初始振型
for i=1:3 %计算三阶固有频率与主振型
pp0=0;
i
B=D*A;
pp=1、0/B(3); %B(3)为B中得最后一个元素
A=B/B(3);
while abs((pp-pp0)/pp)>1、e-6
pp0=pp;
B=D*A;
pp=1、0/B(3);
A=B/B(3);
end
f=sqrt(pp)/2/pi %固有频率单位转换为Hz
fprintf(fid1,'%20、5f',A); %输入主振型数据
fprintf(fid2,'%20、5f',f); %输入固有频率数据
D=D-A*A'*M/(A'*M*A*pp);
end
fid1=fopen('A-1','rt'); %打开主振型文件
A=fscanf(fid1,'%f',[3,3]) %主振型写成矩阵
fid2=fopen('B-1','rt'); %打开固有频率文件
f=fscanf(fid2,'%f',[3,1]) %固有频率写成矩阵
t=1:3;
h1=figure('numbertitle','off','name','0','pos',[50 200 420 420]);
bar(t,f(t,1)),xlabel('频率阶级次'),ylabel('Hz'),
title('固有频率'),hold on,grid
h1=figure('numbertitle','off','name','1','pos',[50 200 420 420]);
bar(t,A(t,1)),xlabel('自由度(质量块)'),ylabel('振型向量'),
title('1阶主振型'),hold on,grid
pause(0、1)
h1=figure('numbertitle','off','name','2','pos',[50 200 420 420]);
bar(t,A(t,2)),xlabel('自由度(质量块)'),ylabel('振型向量'),
title('2阶主振型'),hold on,grid
pause(0、1)
h1=figure('numbertitle','off','name','3','pos',[50 200 420 420]);
bar(t,A(t,3)),xlabel('自由度(质量块)'),ylabel('振型向量'),
title('3阶主振型'),hold on,grid
pause(0、1)
end
%chuandijuzhen、m; %传递矩阵方法求固有频率
clear all,clear close
J1=1;J2=1;J3=2;
k2=1100000;k3=1200000;k4=100000;
fid=fopen('chuandi','wt'); %建立(打开)速度文件
M1L=0;
for WN=0:0、01:2000
shita1R=1;M1R=-WN^2*J1;
shita2R=shita1R+1/k2*M1R;M2R=shita1R*(-WN^2*J2)+(1+(-WN^2*J2)/k2)*M1R;
shita3R=shita2R+1/k3*M2R;M3R=shita2R*(-WN^2*J3)+(1+(-WN^2*J3)/k3)*M2R;
shita4R=shita3R+1/k4*M3R;
if abs(shita4R)<0、005
WN %搜索到得固有频率(rad/s),并显示
shita=[shita1R;shita2R;shita3R;shita4R];%搜索到振型,并显示
bar(shita),xlabel('对应得质量块'),ylabel('振型向量')
pause(1、0)
end
fprintf(fid,'%30、15f',shita4R);
end
fid=fopen('chuandi','rt');
x=fscanf(fid,'%f',[1,200001]);
t=1:200001;
plot(0、01*t,x);grid,xlabel('频率rad/s'),ylabel('第四个质量块得转角(rad/s)'),title('用传递矩阵法求固有频率')
function cdjz2
%chuandijuzhen、m; %传递矩阵方法求固有频率
clear all,clear close
J1=0、5;J2=1;
k2=100000;k3=100000;
fid=fopen('chuandi3','wt'); %建立(打开)速度文件
M1L=0;
for WN=0:0、01:2000
shita1R=1;M1R=-WN^2*J1;
shita2R=shita1R+1/k2*M1R;M2R=shita1R*(-WN^2*J2)+(1+(-WN^2*J2)/k2)*M1R;
shita3R=shita2R+1/k3*M2R;
if abs(shita3R)<0、005
WN %搜索到得固有频率(rad/s),并显示
shita=[shita1R;shita2R;shita3R;] %搜索到振型,并显示
bar(shita),xlabel('对应得质量块'),ylabel('振型向量')
pause(1、0)
end
fprintf(fid,'%30、15f',shita3R);
end
fid=fopen('chuandi3','rt');
x=fscanf(fid,'%f',[1,200001]);
t=1:200001;
plot(0、01*t,x);grid,xlabel('频率rad/s'),ylabel('第四个质量块得转角(rad/s)'),title('用传递矩阵法求固有频率')
function zuoye8
%矩阵迭代法求系统得三阶固有频率与主阵型
clear all
close all
fid1=fopen('A-2','wt'); %建立主振型文件
fid2=fopen('B-2','wt'); %建立固有频率文件
%输入质量矩阵
M(1,1)=1;
M(2,2)=1;
M(3,3)=1;
%输入刚度矩阵
K(1,1)=2;K(1,2)=-1;K(2,1)=- 1;K(2,2)=2;K(2,3)=- 1;K(3,2)=- 1;K(3,3)= 1
%计算特征值与特征向量
D=inv(K)*M; %原始动力矩阵
A=ones(3,1); %初始振型
for i=1:3 %计算三阶固有频率与主振型
pp0=0;
i
B=D*A;
pp=1、0/B(3); %B(3)为B中得最后一个元素
A=B/B(3);
while abs((pp-pp0)/pp)>1、e-6
pp0=pp;
B=D*A;
pp=1、0/B(3);
A=B/B(3);
end
f=sqrt(pp)
fprintf(fid1,'%20、5f',A); %输入主振型数据
fprintf(fid2,'%20、5f',f); %输入固有频率数据
D=D-A*A'*M/(A'*M*A*pp);
end
fid1=fopen('A-2','rt'); %打开主振型文件
A=fscanf(fid1,'%f',[3,3]) %主振型写成矩阵
fid2=fopen('B-2','rt'); %打开固有频率文件
f=fscanf(fid2,'%f',[3,1]) %固有频率写成矩阵
t=1:3;
h1=figure('numbertitle','off','name','0','pos',[50 200 420 420]);
bar(t,f(t,1)),xlabel('频率阶级次'),ylabel('Hz'),
title('固有频率'),hold on,grid
h1=figure('numbertitle','off','name','1','pos',[50 200 420 420]);
bar(t,A(t,1)),xlabel('自由度(质量块)'),ylabel('振型向量'),
title('1阶主振型'),hold on,grid
pause(0、1)
h1=figure('numbertitle','off','name','2','pos',[50 200 420 420]);
bar(t,A(t,2)),xlabel('自由度(质量块)'),ylabel('振型向量'),
title('2阶主振型'),hold on,grid
pause(0、1)
h1=figure('numbertitle','off','name','3','pos',[50 200 420 420]);
bar(t,A(t,3)),xlabel('自由度(质量块)'),ylabel('振型向量'),
title('3阶主振型'),hold on,grid
pause(0、1)
end
function zuoye9
%矩阵迭代法求系统得三阶固有频率与主阵型
clear all
close all
fid1=fopen('A-3','wt'); %建立主振型文件
fid2=fopen('B-3','wt'); %建立固有频率文件
%输入质量矩阵
M(1,1)=4;
M(2,2)=2;
M(3,3)=1;
%输入刚度矩阵
K(1,1)=4;K(1,2)=-1;K(2,1)=- 1;K(2,2)=2;K(2,3)=- 1;K(3,2)=- 1;K(3,3)= 1
%计算特征值与特征向量
D=inv(K)*M; %原始动力矩阵
U=inv(K)
A=ones(3,1); %初始振型
for i=1:3 %计算三阶固有频率与主振型
pp0=0;
i
B=D*A;
pp=1、0/B(1); %B(1)为B中得第一个元素
A=B/B(1);
while abs((pp-pp0)/pp)>1、e-6
pp0=pp;
B=D*A;
pp=1、0/B(1);
A=B/B(1);
end
f=sqrt(pp)
fprintf(fid1,'%20、5f',A); %输入主振型数据
fprintf(fid2,'%20、5f',f); %输入固有频率数据
D=D-A*A'*M/(A'*M*A*pp);
end
fid1=fopen('A-3','rt'); %打开主振型文件
A=fscanf(fid1,'%f',[3,3]) %主振型写成矩阵
fid2=fopen('B-3','rt'); %打开固有频率文件
f=fscanf(fid2,'%f',[3,1]) %固有频率写成矩阵
t=1:3;
h1=figure('numbertitle','off','name','0','pos',[50 200 420 420]);
bar(t,f(t,1)),xlabel('频率阶级次'),ylabel('Hz'),
title('固有频率'),hold on,grid
h1=figure('numbertitle','off','name','1','pos',[50 200 420 420]);
bar(t,A(t,1)),xlabel('自由度(质量块)'),ylabel('振型向量'),
title('1阶主振型'),hold on,grid
pause(0、1)
h1=figure('numbertitle','off','name','2','pos',[50 200 420 420]);
bar(t,A(t,2)),xlabel('自由度(质量块)'),ylabel('振型向量'),
title('2阶主振型'),hold on,grid
pause(0、1)
h1=figure('numbertitle','off','name','3','pos',[50 200 420 420]);
bar(t,A(t,3)),xlabel('自由度(质量块)'),ylabel('振型向量'),
title('3阶主振型'),hold on,grid
pause(0、1)
end
function cdjz2
%chuandijuzhen、m; %传递矩阵方法求固有频率
clear all,clear close
J1=0、5;J2=1;
k2=10000;k3=10000;
fid=fopen('chuandi3','wt'); %建立(打开)速度文件
M1L=0;
for WN=0:0、01:2000
shita1R=1;M1R=-WN^2*J1;
shita2R=shita1R+1/k2*M1R;M2R=shita1R*(-WN^2*J2)+(1+(-WN^2*J2)/k2)*M1R;
shita3R=shita2R+1/k3*M2R;
if abs(shita3R)<0、005
WN %搜索到得固有频率(rad/s),并显示
shita=[shita1R;shita2R;shita3R;] %搜索到振型,并显示
bar(shita),xlabel('对应得质量块'),ylabel('振型向量')
pause(1、0)
end
fprintf(fid,'%30、15f',shita3R);
end
fid=fopen('chuandi3','rt');
x=fscanf(fid,'%f',[1,200001]);
t=1:200001;
plot(0、01*t,x);grid,xlabel('频率rad/s'),ylabel('第三个质量块得转角(rad/s)'),title('用传递矩阵法求固有频率')
第3章
function vtb3(m,c,k,x0,v0,tf,w,f0,delt)
%用欧拉法计算单自由度系统谐迫振动响应
wn=sqrt(k/m); %计算固有频率
fid1=fopen('disp','wt'); %建立一个位移文件disp、dat
for t=0:delt:tf; %delt为时间步长
xdd=(f0*sin(w*t)-k*x0-c*v0)/m; %计算加速度
xd=v0+xdd*delt; %计算速度
x=x0+xd*delt; %计算位移x
fprintf(fid1,'%10、4f',x0); %向文件中写数据
x0=x;v0=xd;
t
end
fid2=fopen('disp','rt'); %打开disp、dat文件
n=tf/delt; %disp、dat文件中位移得个数
x=fscanf(fid2,'%f',[1,n]); %将disp、dat文件中文艺写成矩阵
t=1:n;
plot(t,x),grid
xlabel('时间(s)')
ylabel('位移(s)')
title('位移与时间得关系')
function vtb4(m,c,k,x0,v0,tf,w,f0,delt)
%用改进得欧拉法计算单自由度系统谐迫振动响应
wn=sqrt(k/m); %计算固有频率
fid1=fopen('disp','wt'); %建立一个位移文件disp、dat
for t=0:delt:tf; %delt为时间步长
xdd=(f0*sin(w*t)-k*x0-c*v0)/m; %计算加速度
x3d=(f0*w*cos(w*t)-k*v0-c*xdd)/m;
xd=v0+xdd*delt+x3d*delt^2/2; %计算速度
x=x0+xd*delt+xdd*delt^2/2; %计算位移x
fprintf(fid1,'%10、4f',x0); %向文件中写数据
x0=x;v0=xd;
t
end
fid2=fopen('disp','rt'); %打开disp、dat文件
n=tf/delt; %disp、dat文件中位移得个数
x=fscanf(fid2,'%f',[1,n]); %将disp、dat文件中文艺写成矩阵
t=1:n;
plot(t,x),grid
xlabel('时间(s)')
ylabel('位移(s)')
title('位移与时间得关系')
function vtb5(tf,delt) %用线性加速度法计算三自由度系统谐迫振动响应,tf为仿真时间,delt为仿真时间步长delt
close all;clc
fid1=fopen('disp5','wt'); %建立一个位移文件dip5、dat
m=2*[1 0 0;0 1 0;0 0 1]; %质量矩阵
c=1、5*[2 -1 0;-1 2 -1;0 -1 2]; %阻尼矩阵
k=50*[2 -1 0;-1 2 -1;0 -1 2]; %刚度矩阵
x0=[1 1 1]'; %初始位移
v0=[1 1 1]'; %初始速度
md=inv(m+delt/2*c+1/6*delt^2*k);
m1=inv(m);
[E,F]=eig(m1*k);
diag(sqrt(F)); %计算固有频率(rad/s)
for t=0:delt:tf; %delt为时间步长
f=[2、00*sin(3、754*t) -2、00*cos(2、2*t) 1、00*sin(2、8*t)]';
if t==0; xdd0=m1*(f-k*x0-c*v0); %计算初始加速度
else
x=md*(m*(x0+delt*v0+delt^2/3*xdd0)+c*(delt/2*x0+delt^2/3*v0+delt^3/12*xdd0)+delt^2/6*f);%计算位移
xdd=6/delt^2*(x-(x0+delt*v0+delt^2/3*xdd0)); %计算加速度
xd=v0+delt/2*(xdd0+xdd);%计算速度
xdd0=xdd;v0=xd;x0=x;
fprintf(fid1,'%10、4f',x0);%向文件中写数据
t %显示计算时间步长
end
end
fid2=fopen('disp5','rt'); %打开disp5、dat文件
n=tf/delt; %disp5、dat文件中位移得个数
x=fscanf(fid2,'%f',[3,n]); %将disp5、dat文件中得位移写成矩阵
t=1:n;
figure('numbertitle','off','name','自由度-1得位移','pos',[450 180 400 420]);
plot(t,x(1,t)),grid,xlabel('时间*0、1秒'),title('自由度-1得位移与时间得关系')
figure('numbertitle','off','name','自由度-2得位移','pos',[350 160 400 420]);
plot(t,x(2,t)),grid,xlabel('时间*0、1秒'),title('自由度-3得位移与时间得关系')
figure('numbertitle','off','name','自由度-3得位移','pos',[250 140 400 420]);
plot(t,x(3,t)),grid,xlabel('时间*0、1秒'),title('自由度-3得位移与时间得关系')
function vtb6(tf,delt) %用线纽马克-β法计算三自由度系统谐迫振动响应,tf为仿真时间,delt为仿真时间步长delt
close all;clc
fid1=fopen('disp6','wt'); %建立一个位移文件dip6、dat
m=2*[1 0 0;0 1 0;0 0 1]; %质量矩阵
c=1、5*[2 -1 0;-1 2 -1;0 -1 2]; %阻尼矩阵
k=50*[2 -1 0;-1 2 -1;0 -1 2]; %刚度矩阵
x0=[1 1 1]'; %初始位移
v0=[1 1 1]'; %初始速度
bita=1/6;
md=inv(m+delt/2*c+bita*delt^2*k);
m1=inv(m);
[E,F]=eig(m1*k);
diag(sqrt(F)); %计算固有频率(rad/s)
for t=0:delt:tf; %delt为时间步长
f=[2、00*sin(3、754*t) -2、00*cos(2、2*t) 1、00*sin(2、8*t)]';
if t==0; xdd0=m1*(f-k*x0-c*v0); %计算初始加速度
else
xdd=md*(f-c*(v0+delt/2*xdd0)-k*(x0+delt*v0+(1/2-bita)*delt^2*xdd0)); %计算加速度
xd=v0+delt/2*(xdd0+xdd);%计算速度
x=x0+delt*v0+delt^2/2*xdd0+bita*delt^3*(xdd-xdd0)/delt; %计算位移
v0=xd;x0=x;
fprintf(fid1,'%10、4f',x0);%向文件中写数据
t %显示计算时间步长
end
end
fid2=fopen('disp6','rt'); %打开disp6、dat文件
n=tf/delt; %disp6、dat文件中位移得个数
x=fscanf(fid2,'%f',[3,n]); %将disp6、dat文件中得位移写成矩阵
t=1:n;
figure('numbertitle','off','name','自由度-1得位移','pos',[450 180 400 420]);
plot(t,x(1,t)),grid,xlabel('时间*0、1秒'),title('自由度-1得位移与时间得关系')
figure('numbertitle','off','name','自由度-2得位移','pos',[350 160 400 420]);
plot(t,x(2,t)),grid,xlabel('时间*0、1秒'),title('自由度-2得位移与时间得关系')
figure('numbertitle','off','name','自由度-3得位移','pos',[250 140 400 420]);
plot(t,x(3,t)),grid,xlabel('时间*0、1秒'),title('自由度-3得位移与时间得关系')
function vtb7(tf,delt) %用威尔逊θ法计算三自由度系统谐迫振动响应,tf为仿真时间,delt为仿真时间步长delt
close all;clc
fid1=fopen('disp7','wt'); %建立一个位移文件dip6、dat
m=2*[1 0 0;0 1 0;0 0 1]; %质量矩阵
c=1、5*[2 -1 0;-1 2 -1;0 -1 2]; %阻尼矩阵
k=50*[2 -1 0;-1 2 -1;0 -1 2]; %刚度矩阵
x0=[1 1 1]'; %初始位移
v0=[1 1 1]'; %初始速度
theta=1、4;
md=inv(k+3*c/theta/delt+6/(theta*delt^2)*m);
m1=inv(m);
[E,F]=eig(m1*k);
diag(sqrt(F)); %计算固有频率(rad/s)
for t=0:delt:tf; %delt为时间步长
f=[2、00*sin(3、754*t) -2、00*cos(2、2*t) 1、00*sin(2、8*t)]';
if t==0; xdd0=m1*(f-k*x0-c*v0); %计算初始加速度
else
xtheta=md*(m*(2*xdd0+6/theta/delt*v0+6/(theta*delt)^2*x0)+c*(theta*delt/2*xdd0+2*v0+3/theta/delt*x0)+f); %计算(t+θdelt)时刻得速度
xddtheta=6/(theta*delt)^2*(xtheta-x0)-6/theta/delt*v0-2*xdd0; %计算(t+θdelt)时刻得加速度
xdd=(1-1/theta)*xdd0+1/theta*xddtheta; %计算(t+delt)时刻得加速度
xd=v0+delt/2*(xdd0+xdd); %计算(t+delt)速度
x=x0+delt*v0+delt^2*(2*xdd0+xdd)/6; %计算(t+delt)位移
v0=xd;x0=x;xdd0=xdd;
fprintf(fid1,'%10、4f',x0);%向文件中写数据
t %显示计算时间步长
end
end
fid2=fopen('disp7','rt'); %打开disp6、dat文件
n=tf/delt; %disp7、dat文件中位移得个数
x=fscanf(fid2,'%f',[3,n]); %将disp7、dat文件中得位移写成矩阵
t=1:n;
figure('numbertitle','off','name','自由度-1得位移','pos',[450 180 400 420]);
plot(t,x(1,t)),grid,xlabel('时间*0、1秒'),title('自由度-1得位移与时间得关系')
figure('numbertitle','off','name','自由度-2得位移','pos',[350 160 400 420]);
plot(t,x(2,t)),grid,xlabel('时间*0、1秒'),title('自由度-2得位移与时间得关系')
figure('numbertitle','off','name','自由度-3得位移','pos',[250 140 400 420]);
plot(t,x(3,t)),grid,xlabel('时间*0、1秒'),title('自由度-3得位移与时间得关系')