《动力学第三章(13页).doc》由会员分享，可在线阅读，更多相关《动力学第三章(13页).doc（13页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、-第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-
2、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
3、(固有频率为%.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=fsc
4、anf(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(振型向量
5、),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)e
6、nd%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
7、+(-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);
8、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*
9、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(
10、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,
11、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=fop
12、en(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
13、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(n
14、umbertitle,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)
15、,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=f
16、open(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(振型向量) p
17、ause(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
18、.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)yla
19、bel(位移(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
20、/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
21、;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
22、.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 %显示计算时
23、间步长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
24、(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
25、 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
26、*(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; %dis
27、p6.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(n
28、umbertitle,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;-
29、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
30、/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
31、)/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的位移与时间的关系)-第 13 页-