动力学程序.doc

动力学程序 第一章 %第8页，有阻尼单自由度问题 function VTB1(m,c,k,x0,v0,tf) %m是质量，c,k刚度，x0质量块位移，v0质量块速度 clc %清屏 wn=sqrt(k/m); z=c/2/m/wn;%ζ=n/wn,n=c/2m wd=wn*sqrt(1-z^2); fprintf(固有频率为%.3g.rad/s.\n,wn);%输出wn fprintf(阻尼系数为%.3g.\n,z);%输出ξ fprintf(有阻尼的固有频率为%.3g.rad/s.\n,wd);%输出wd t=0:tf/10000:tf;%时域的长短，决定很坐标的长短 if z1 %判断ζ的取值，如果ζ1,弱阻尼情况，按以下1-10到1-12公式 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);%输出A fprintf(phi=%.3g\n,phi);%输出ψ elseif z==1 %临界阻尼情况，按照1-14到1-15公式计算 a1=x0;% C1 a2=v0+wn*x0;%C2 fprintf(a1=%.3g\n,a1);%输出a1 fprintf(a2=%.3g\n,a2);%输出a2 x=(a1+a2*t).*exp(-wn*t); else %过阻尼，按照1-13公式计算 a1=(-v0+(-z+sqrt(z^2-1))*wn*x0)/2/wn/sqrt(z^2-1); a2=(v0+(z+sqrt(z^2-1))*wn*v0)/2/wn/sqrt(z^2-1); fprintf(a1=%.3g\n,a1);%输出a1 fprintf(a2=%.3g\n,a2);%输出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 on xlabel(时间(s)) ylabel(位移(m)) title(位移相对时间的关系) 单自由度系统的谐迫振动（P11业 function vtb2(m,c,k,x0,v0,tf,w,f0) %单自由度系统的谐迫振动 clc %清屏 wn=sqrt(k/m); z=c/2/m/wn;%ζ=n/wn,n=c/2m lan=w/wn %λ的求法1-18公式 wd=wn*sqrt(1-z^2); fprintf(固有频率为%.3g.rad/s.\n,wn); %输出Wn fprintf(阻尼系数为%.3g.\n,z);%输出ζ fprintf(有阻尼的固有频率为%.3g.rad/s.\n,wd); %输出Wd a=sqrt(((v0+z*wn*x0)^2+(x0*wd)^2)/wd^2); t=0:tf/1000:tf; phi1=atan(x0*wd/(v0+z*wn*x0));%按1-12求ψ phi2=atan(2*z*lan/(1-lan^2));%求Φ b=wn^2*f0/k/sqrt((wn^2-w^2)+(2*z*wn*w)^2);%求B的稳态响应的振幅 x=a*exp(-z*wn*t).*sin(sqrt(1-z^2)*wn*t+phi1)+b*sin(w*t-phi2);%响应方程 plot(real(t),real(x)),grid xlabel(时间(s)) ylabel(位移(cm)) title(位移与时间的关系) 第二章 一矩阵迭代法38业例3） function jzdd %矩阵迭代法 clear all %清空当前所有数据 close all %关闭当前所有的绘图窗口 fid1=fopen(A-1,wt); %建立第一个名为“A-1”的可写文档 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; %inv表示对矩阵取逆，公式2-65 A=ones(3,1);%定义一个初始迭代阵型，ones()函数表示是3个位为1的单列矩阵，ones(i,j) %则是i行两列都是j都是1的数组！在这方法中一般取ones(i,1),i=质量各数 fo