動力學第三章

1、第2章function VTB2(m,c,k,x0,v0,tf,w,f0)%單自由度系統(tǒng)的諧迫振動clcwn=sqrt(k/m);z=c/2/m/wn;lan=w/wnwd=w n*sqrt(1-zA2);A=sqrt(vO+z*w n*xO)A2+(xO*wdf2)/wdA2);t=0:tf/1000:tf;phi1=atan2(x0*wd,v0+z*wn*x0);phi2=ata n2(2*z*la n,1-la 門人2);B=w nA2*fO/k/sqrt(w nA2-wA2F2+(2*z*w n*wF2);x=A*exp(-z*w n*t).*si n( sqrt(1-zA2)*w n

2、*t+phi1)+B*si n(w*t-phi2);plot(t,x),gridxlabel( 時間 (s) )ylabel( 位移 )title( 位移與時間的關系 )function VTB1(m,c,k,x0,v0,tf)%VTB傭來計算單自由度有阻尼自由振動系統(tǒng)的響應%VTB1繪岀單自由度有阻尼自由振動系統(tǒng)的響應圖%n為質量;c為阻尼;k為剛度;x0為初始位移;v0為初始速度;tf為仿真時間%VTB1(zeta,w,x0,v0,tf) 繪出單自由度有阻尼自由振動系統(tǒng)的響應圖%程序中z為阻尼系數(shù)E ;wn為固有頻率3 n;A為振動幅度；phi為初相位0 clcwn=sqrt(k/m);z

3、=c/2/m/wn;wd=w n*sqrt(1-zA2);fprintf( 固有頻率為 %.s.n,wd);fprintf( 阻尼系數(shù)為 %.3g.n,z);fprintf( 有阻尼的固有頻率為%.s.n ,wd);t=0:tf/1000:tf;if zpp0=pp;B=D*A;pp=B(3);A=B/B(3);endf=sqrt(pp)/2/pi%固有頻率單位轉換為Hzfprintf(fid1,% ,A);%輸入主振型數(shù)據(jù)fprintf(fid2,% ,f);%輸入固有頻率數(shù)據(jù)D=D-A*A*M/(A*M*A*pp);endfid1=fopen( A-1, rt );%打開主振型文件A=fs

4、canf(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( 固有頻率 ),holdon,gridh1=figure( numbertitle , off , name ,1 , pos ,50 200 420 420);bar(t,A(t,1),xla

5、bel( 自由度(質量塊) ),ylabel( 振型向量 ), title( 1 階主振型 ),hold on,gridpauseh1=figure( numbertitle , off , name ,2 , pos ,50 200 420 420); bar(t,A(t,2),xlabel( 自由度(質量塊) ),ylabel( 振型向量 ), title( 2 階主振型 ),hold on,gridpauseh1=figure( numbertitle , off , name ,3 , pos ,50 200 420 420); bar(t,A(t,3),xlabel( 自由度(質量塊

6、) ),ylabel( 振型向量 ), title( 3 階主振型 ),hold on,gridpauseend%; %傳遞矩陣方法求固有頻率clear all ,clear closeJ1=1;J2=1;J3=2;k2=1100000;k3=1200000;k4=100000;fid=fopen( chuandi , wt ); %建立(打開)速度文件M1L=0;for WN=0:2000shita1R=1;M1R=-WNA2*J1;shita2R=shita1R+1/k2*M1R;M2R=shita1R*(-WNA2*J2)+(1+(-WNA2*J2)/k2)*M1R;shita3R=sh

7、ita2R+1/k3*M2R;M3R=shita2R*(-WNA2*J3)+(1+(-WNA2*J3)/k3)*M2R;shita4R=shita3R+1/k4*M3R;if abs(shita4R)WN%搜索到的固有頻率 (rad/s), 并顯示shita=shita1R;shita2R;shita3R;shita4R;%搜索到振型 , 并顯示bar(shita),xlabel( 對應的質量塊 ),ylabel( 振型向量 ) pauseendfprintf(fid, % ,shita4R);endfid=fopen( chuandi , rt );x=fscanf(fid, %f ,1,2

8、00001);t=1:200001;plot*t,x);grid,xlabel( 頻率 rad/s ),ylabel( 第四個質量塊的轉角 (rad/s) ),title( 用傳遞矩 陣法求固有頻率 )function cdjz2%; %傳遞矩陣方法求固有頻率clear all ,clear closeJ1=;J2=1;k2=100000;k3=100000;fid=fopen( chuandi3 , wt ); %建立(打開)速度文件M1L=0;for WN=0:2000shita1R=1;M1R=-WNA2*J1;shita2R=shita1R+1/k2*M1R;M2R=shita1R*(

9、-WNA2*J2)+(1+(-WNA2*J2)/k2)*M1R;shita3R=shita2R+1/k3*M2R;if abs(shita3R)pp0=pp;B=D*A;pp=B(3);A=B/B(3);endf=sqrt(pp)fprintf(fid1,% ,A);%輸入主振型數(shù)據(jù)fprintf(fid2,% ,f);%輸入固有頻率數(shù)據(jù)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(fid

10、2, %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( 固有頻率 ),holdon,gridh1=figure( numbertitle, off , name ,1, pos,50 200 420 420);bar(t,A(t,1),xlabel( 自由度（質量塊） ),ylabel( 振型向量 ）,title( 1 階主振型 ),holdon,gridpauseh1=figure(

11、numbertitle, off , name ,2, pos,50 200 420 420);bar(t,A(t,2),xlabel( 自由度（質量塊） ),ylabel( 振型向量 ）,title( 2 階主振型 ),holdon,gridpauseh1=figure( numbertitle, off , name ,3, pos,50 200 420 420);bar(t,A(t,3),xlabel( 自由度（質量塊） ),ylabel( 振型向量 ）,title( 3 階主振型 ),holdon,gridpauseendfunction zuoye9%矩陣迭代法求系統(tǒng)的三階固有頻率和

12、主陣型 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;iB=D*A;pp=B(1);%B(1)為B

13、中的第一個元素A=B/B(1);while abs(pp-pp0)/pp)pp0=pp;B=D*A;pp=B(1);A=B/B(1);endf=sqrt(pp)fprintf(fid1, % ,A); %輸入主振型數(shù)據(jù)fprintf(fid2, % ,f);%輸入固有頻率數(shù)據(jù)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=fig

14、ure( numbertitle, off , name ,0 , pos,50 200 420 420);bar(t,f(t,1),xlabel( 頻率階級次 ),ylabel( Hz),title( 固有頻率 ),holdon,gridh1=figure( numbertitle, off , name ,1 , pos,50 200 420 420);bar(t,A(t,1),xlabel( 自由度(質量塊) ),ylabel( 振型向量 ),title( 1 階主振型 ),holdon,gridpauseh1=figure( numbertitle, off , name ,2 , p

15、os,50 200 420 420);bar(t,A(t,2),xlabel( 自由度(質量塊) ),ylabel( 振型向量 ),title( 2 階主振型 ),holdon,gridpauseh1=figure( numbertitle, off , name ,3 , pos,50 200 420 420);bar(t,A(t,3),xlabel( 自由度(質量塊) ),ylabel( 振型向量 ),title( 3 階主振型 ),holdon,gridpauseend function cdjz2%; %傳遞矩陣方法求固有頻率clear all ,clear closeJ1=;J2=1

16、;k2=10000;k3=10000;fid=fopen( chuandi3 , wt ); %建立(打開)速度文件M1L=0;for WN=0:2000shita1R=1;M1R=-WNA2*J1; shita2R=shita1R+1/k2*M1R;M2R=shita1R*(-WNA2*J2)+(1+(-WNA2*J2)/k2)*M1R;shita3R=shita2R+1/k3*M2R;if abs(shita3R)WN%搜索到的固有頻率 (rad/s), 并顯示shita=shita1R;shita2R;shita3R;%搜索到振型 , 并顯示bar(shita),xlabel(對應的質量

17、塊 ),ylabel( 振型向量 )pauseendfprintf(fid, % ,shita3R);endfid=fopen( chuandi3 , rt );x=fscanf(fid, %f ,1,200001);t=1:200001;plot*t,x);grid,xlabel(頻率 rad/s ),ylabel( 第三個質量塊的轉角 (rad/s) ),title( 用傳遞矩陣法求固有頻率 )第3章function vtb3(m,c,k,x0,v0,tf,w,f0,delt)%用歐拉法計算單自由度系統(tǒng)諧迫振動響應wn=sqrt(k/m);%計算固有頻率fid1=fopen(disp ,

18、wt );%建立一個位移文件for t=0:delt:tf;%delt 為時間步長xdd=(f0*sin(w*t)-k*x0-c*v0)/m;%計算加速度xd=v0+xdd*delt;%計算速度x=x0+xd*delt;%計算位移 xfprintf(fid1, % ,x0);%向文件中寫數(shù)據(jù)x0=x;v0=xd;tendfid2=fopen( disp , rt );%打開文件n=tf/delt;%文件中位移的個數(shù)t=1:n;plot(t,x),gridxlabel( 時間 (s) )ylabel( 位移 (s) )title( 位移與時間的關系 )function vtb4(m,c,k,x0

19、,v0,tf,w,f0,delt)%用改進的歐拉法計算單自由度系統(tǒng)諧迫振動響應 wn=sqrt(k/m);%計算固有頻率fid1=fopen( disp , wt );%建立一個位移文件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=vO+xdd*delt+x3d*deM2/2;算速度x=xO+xd*delt+xdd*deM2/2;算位移 xfprintf(fid1,% ,x0);%向文件中寫數(shù)據(jù)x0=x;v0=xd;tendfid2=fop

20、en( disp , rt );%打開文件n=tf/delt;%文件中位移的個數(shù)x=fscanf(fid2, %f ,1,n);%將文件中文藝寫成矩陣t=1:n;plot(t,x),gridxlabel( 時間 (s) )ylabel( 位移 (s) )title( 位移與時間的關系 )function vtb5(tf,delt)%用線性加速度法計算三自由度系統(tǒng)諧迫振動響應，tf 為仿真時間 ,delt為仿真時間步長 deltclose all ;clcfid1=fopen( disp5 , wt ); %建立一個位移文件m=2*1 0 0;0 1 0;0 0 1;%質量矩陣c=*2 -1 0

21、;-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=i nv(m+delt/2*c+1/6*deltA2*k);m1=inv(m);E,F=eig(m1*k);diag(sqrt(F);%計算固有頻率( rad/s )for t=0:delt:tf; %delt 為時間步長f=*sin*t) *cos*t) *sin*t);if t=0; xdd0=m1*(f-k*x0-c*v0);%計算初始加速度elsex=md*(m*(x0+delt*v0+deltA2/3*xdd0)+c*

22、(delt/2*x0+deltA2/3*v0+deltA3/12*xdd0)+deltA2/6*f); %計算位移xdd=6/deltA2*(x-(x0+delt*v0+deltA2/3*xdd0);%計算加速度xd=v0+delt/2*(xdd0+xdd);%計算速度xdd0=xdd;v0=xd;x0=x;fprintf(fid1, % ,x0); %向文件中寫數(shù)據(jù) t %顯示計算時間步長endendfid2=fopen( disp5 , rt ); %打開文件 n=tf/delt;%文件中位移的個數(shù)x=fscanf(fid2,%f ,3,n);%將文件中的位移寫成矩陣t=1:n;figur

23、e( numbertitle, off,name , 自由度 -1的位移 , pos ,450 180 400 420);plot(t,x(1,t),grid,xlabel( 時間 *秒 ),title(自由度 -1 的位移與時間的關系 )figure( numbertitle, off,name , 自由度 -2的位移 , pos ,350 160 400 420);plot(t,x(2,t),grid,xlabel( 時間 *秒 ),title(自由度 -3 的位移與時間的關系 )figure( numbertitle, off,name , 自由度 -3的位移 , pos ,250 14

24、0 400 420);plot(t,x(3,t),grid,xlabel( 時間 *秒 ),title(自由度 -3 的位移與時間的關系 )function vtb6(tf,delt)%用線紐馬克-B法計算三自由度系統(tǒng)諧迫振動響應，tf為仿真時間，delt為仿真時間步長 delt close all ;clcfid1=fopen( disp6 , wt );%建立一個位移文件m=2*1 0 0;0 1 0;0 0 1;%質量矩陣c=*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

25、;%初始速度 bita=1/6;md=i nv(m+delt/2*c+bita*deltA2*k);m1=inv(m);E,F=eig(m1*k);diag(sqrt(F);%計算固有頻率( rad/s )for t=0:delt:tf; %delt 為時間步長f=*sin*t) *cos*t) *sin*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)*deltA2*xdd0);%計算加速度xd=v0+delt/2*(xdd0+xdd);%計算速

26、度x=x0+delt*v0+deltA2/2*xdd0+bita*deltA3*(xdd-xdd0)/delt;%計算位移v0=xd;x0=x;fprintf(fid1, % ,x0); %向文件中寫數(shù)據(jù)t%顯示計算時間步長endendfid2=fopen( disp6 , rt );%打開文件n=tf/delt;%文件中位移的個數(shù)x=fscanf(fid2,%f ,3,n);%將文件中的位移寫成矩陣t=1:n;figure( numbertitle, off,name ,自由度 -1 的位移 , pos ,450 180 400 420);plot(t,x(1,t),grid,xlabel(

27、 時間 *秒 ）,title（ 自由度 -1 的位移與時間的關系 ）figure( numbertitle , off ,name ,自由度 -2 的位移 , pos ,350 160 400 420);plot(t,x(2,t),grid,xlabel( 時間 *秒 ),title( 自由度 -2 的位移與時間的關系 )figure( numbertitle , off , name , 自由度 -3 的位移 , pos ,250 140 400 420);plot(t,x(3,t),grid,xlabel( 時間 *秒 ),title( 自由度 -3 的位移與時間的關系 )function

28、 vtb7(tf,delt)%用威爾遜0法計算三自由度系統(tǒng)諧迫振動響應,tf 為仿真時間 ,delt 為仿真時間步長 deltclose all ;clcfid1=fopen( disp7 , wt); %建立一個位移文件m=2*1 0 0;0 1 0;0 0 1;%質量矩陣c=*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=;md=i nv(k+3*c/theta/delt+6/(theta*deltA2)*m);m1=inv(m);E,F=eig(m1*k);diag(sqrt(F);%計算固有頻率( rad/s )for t