动力学第三章.doc

第2章function VTB2(m,c,k,x0,v0,tf,w,f0)%单自由度系统得谐迫振动clcwn=sqrt(k/m);z=c/2/m/wn;lan=w/wnwd=wn*sqrt(1-z2);A=sqrt(v0+z*wn*x0)2+(x0*wd)2)/wd2);t=0:tf/1000:tf;phi1=atan2(x0*wd,v0+z*wn*x0);phi2=atan2(2*z*lan,1-lan2);B=wn2*f0/k/sqrt(wn2-w2)2+(2*z*wn*w)2);x=A*exp(-z*wn*t)、*sin(sqrt(1-z2)*wn*t+phi1)+B*sin(w*t-phi2);plot(t,x),gridxlabel(时间(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为初相位clcwn=sqrt(k/m);z=c/2/m/wn;wd=wn*sqrt(1-z2);fprintf(固有频率为%、3g、rad/s、n,wd);fprintf(阻尼系数为%、3g、n,z);fprintf(有阻尼得固有频率为%、3g、rad/s、n,wd);t=0:tf/1000:tf;if z1、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);endfid1=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,gridh1=figure(numbertitle,off,name,1,pos,50 200 420 420);bar(t,A(t,1),xlabel(自由度(质量块),ylabel(振型向量),title(1阶主振型),hold on,gridpause(0、1)h1=figure(numbertitle,off,name,2,pos,50 200 420 420);bar(t,A(t,2),xlabel(自由度(质量块),ylabel(振型向量),title(2阶主振型),hold on,gridpause(0、1)h1=figure(numbertitle,off,name,3,pos,50 200 420 420);bar(t,A(t,3),xlabel(自由度(质量块),ylabel(振型向量),title(3阶主振型),hold on,gridpause(0、1)end%chuandijuzhen、m; %传递矩阵方法求固有频率clear all,clear closeJ1=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=-WN2*J1; shita2R=shita1R+1/k2*M1R;M2R=shita1R*(-WN2*J2)+(1+(-WN2*J2)/k2)*M1R; shita3R=shita2R+1/k3*M2R;M3R=shita2R*(-WN2*J3)+(1+(-WN2*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);endfid=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 closeJ1=0、5;J2=1;k2=100000;k3=100000;fid=fopen(chuandi3,wt); %建立(打开)速度文件M1L=0;for WN=0:0、01:2000 shita1R=1;M1R=-WN2*J1; shita2R=shita1R+1/k2*M1R;M2R=shita1R*(-WN2*J2)+(1+(-WN2*J2)/k2)*M1R; shita3R=shita2R+1/k3*M2R; if abs(shita3R)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);endfid1=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,gridh1=figure(numbertitle,off,name,1,pos,50 200 420 420);bar(t,A(t,1),xlabel(自由度(质量块),ylabel(振型向量),title(1阶主振型),hold on,gridpause(0、1)h1=figure(numbertitle,off,name,2,pos,50 200 420 420);bar(t,A(t,2),xlabel(自由度(质量块),ylabel(振型向量),title(2阶主振型),hold on,gridpause(0、1)h1=figure(numbertitle,off,name,3,pos,50 200 420 420);bar(t,A(t,3),xlabel(自由度(质量块),ylabel(振型向量),title(3阶主振型),hold on,gridpause(0、1)endfunction zuoye9%矩阵迭代法求系统得三阶固有频率与主阵型clear allclose allfid1=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);endfid1=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,gridh1=figure(numbertitle,off,name,1,pos,50 200 420 420);bar(t,A(t,1),xlabel(自由度(质量块),ylabel(振型向量),title(1阶主振型),hold on,gridpause(0、1)h1=figure(numbertitle,off,name,2,pos,50 200 420 420);bar(t,A(t,2),xlabel(自由度(质量块),ylabel(振型向量),title(2阶主振型),hold on,gridpause(0、1)h1=figure(numbertitle,off,name,3,pos,50 200 420 420);bar(t,A(t,3),xlabel(自由度(质量块),ylabel(振型向量),title(3阶主振型),hold on,gridpause(0、1)endfunction cdjz2%chuandijuzhen、m; %传递矩阵方法求固有频率clear all,clear closeJ1=0、5;J2=1;k2=10000;k3=10000;fid=fopen(chuandi3,wt); %建立(打开)速度文件M1L=0;for WN=0:0、01:2000 shita1R=1;M1R=-WN2*J1; shita2R=shita1R+1/k2*M1R;M2R=shita1R*(-WN2*J2)+(1+(-WN2*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);endfid=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、datfor 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;tendfid2=fopen(disp,rt); %打开disp、dat文件n=tf/delt; %disp、dat文件中位移得个数x=fscanf(fid2,%f,1,n); %将disp、dat文件中文艺写成矩阵t=1:n;plot(t,x),gridxlabel(时间(s)ylabel(位移(s)title(位移与时间得关系)function vtb4(m,c,k,x0,v0,tf,w,f0,delt)%用改进得欧拉法计算单自由度系统谐迫振动响应wn=sqrt(k/m); %计算固有频率fid1=fopen(disp,wt); %建立一个位移文件disp、datfor 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*delt2/2; %计算速度 x=x0+xd*delt+xdd*delt2/2; %计算位移x fprintf(fid1,%10、4f,x0); %向文件中写数据x0=x;v0=xd;tendfid2=fopen(disp,rt); %打开disp、dat文件n=tf/delt; %disp、dat文件中位移得个数x=fscanf(fid2,%f,1,n); %将disp、dat文件中文艺写成矩阵t=1:n;plot(t,x),gridxlabel(时间(s)ylabel(位移(s)title(位移与时间得关系)function vtb5(tf,delt) %用线性加速度法计算三自由度系统谐迫振动响应,tf为仿真时间,delt为仿真时间步长deltclose all;clcfid1=fopen(disp5,wt); %建立一个位移文件dip5、datm=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*delt2*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); %计算初始加速度elsex=md*(m*(x0+delt*v0+delt2/3*xdd0)+c*(delt/2*x0+delt2/3*v0+delt3/12*xdd0)+delt2/6*f);%计算位移xdd=6/delt2*(x-(x0+delt*v0+delt2/3*xdd0); %计算加速度xd=v0+delt/2*(xdd0+xdd);%计算速度xdd0=xdd;v0=xd;x0=x;fprintf(fid1,%10、4f,x0);%向文件中写数据t %显示计算时间步长endendfid2=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为仿真时间步长deltclose all;clcfid1=fopen(disp6,wt); %建立一个位移文件dip6、datm=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*delt2*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); %计算初始加速度elsexdd=md*(f-c*(v0+delt/2*xdd0)-k*(x0+delt*v0+(1/2-bita)*delt2*xdd0); %计算加速度xd=v0+delt/2*(xdd0+xdd);%计算速度x=x0+delt*v0+delt2/2*xdd0+bita*delt3*(xdd-xdd0)/delt; %计算位移v0=xd;x0=x;fprintf(fid1,%10、4f,x0);%向文件中写数据t %显示计算时间步长endendfid2=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为仿真时间步长deltclose all;clcfid1=fopen(disp7,wt); %建立一个位移文件dip6、datm=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*delt2)*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); %计算初始加速度elsextheta=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+delt2*(2*xdd0+xdd)/6; %计算(t+delt)位移v0=xd;x0=x;xdd0=xdd;fprintf(fid1,%10、4f,x0);%向文件中写数据t %显示计算时间步长endendfid2=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得位移与时间得关系)