动力学程序讲解

动力学程序（附程序的详细翻译）第一章第8页，有阻尼单自由度问题%第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/2mwd=wn*sqrt(1-z^2);fprintf('固有频率为%.3g.rad/s.\n',wn);%输出wnfprintf('阻尼系数为%.3g.\n',z);%输出ξfprintf('有阻尼的固有频率为%.3g.rad/s.\n',wd);%输出wdt=0:tf/10000:tf;%时域的长短，决定很坐标的长短if z<1 %判断ζ的取值，如果ζ<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);%输出Afprintf('phi=%.3g\n',phi);%输出ψelseif z==1 %临界阻尼情况，按照1-14到1-15公式计算a1=x0;% C1a2=v0+wn*x0;%C2fprintf('a1=%.3g\n',a1);%输出a1fprintf('a2=%.3g\n',a2);%输出a2x=(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);%输出a1fprintf('a2=%.3g\n',a2);%输出a2x=exp(-z*wn*t).*(a1*exp(-wn*sqrt(z^2-1)*t)+a2*exp(wn*sqrt(z^2-1)*t));endplot(t,x),grid onxlabel('时间(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/2mlan=w/wn %λ的求法1-18公式wd=wn*sqrt(1-z^2);fprintf('固有频率为%.3g.rad/s.\n',wn); %输出Wnfprintf('阻尼系数为%.3g.\n',z);%输出ζfprintf('有阻尼的固有频率为%.3g.rad/s.\n',wd); %输出Wda=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)),gridxlabel('时间(s)')ylabel('位移(cm)')title('位移与时间的关系')第二章一、矩阵迭代法（P38业，例2-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-65A=ones(3,1);%定义一个初始迭代阵型，ones()函数表示是3个位为1的单列矩阵，ones(i,j)%则是i行两列都是j都是1的数组！在这方法中一般取ones(i,1),i=质量各数for i=1:3 %进入for循环，将要执行书本上37业的迭代，循环次数为1到i，没执行一次循环代表一个质量块的主阵型pp0=0; %令初始的P（0）=0，并且输出Ｐ（０）i %输出iB=D*A; %代入2-60的公式pp=1.0/B(3); %同上，B(3),是你当时定义主阵型时，第i个元素全为1，当然也可以是第一个A=B/B(3); %同上%以上三段代码，是初始赋值while abs((pp-pp0)/pp)>1e-6 % 判断是否满足迭代的条件|P（k）^2-P(k-1)^2|<δ%满足后退出while循环pp0=pp; %P(k)^2=上诉最终的pp值B=D*A;pp=1.0/B(3);A=B/B(3);%真正地迭代代码end %结束while循环f=sqrt(pp)/2/pi %公式2-61fprintf(fid1,'%20.5f',A); %fid为文件句柄，指定要写入数据的文件，format是用来控制所写数据格式的格式符，与fscanf函数相同，A是用来存放数据的矩阵。

fprintf(fid2,'%20.5f',f); %储存fD=D-A*A'*M/(A'*M*A*pp); %下一阶的动力矩阵【D*】，A’代表A的逆矩阵end%结束for循环fid1=fopen('A-1','rt');%读取A-1文档A=fscanf(fid1,'%f',[3,3]);%输入3x3的A矩阵，来源为fid1fid2=fopen('B-1','rt');%读取B-1文档f=fscanf(fid2,'%f',[3,1]);%输入（拾取）3x1的f矩阵,来源是fid2t=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 %bar代表绘制长条图t代表起点横坐标，f(t,1)纵坐标,输出各阶频率阶次的的固有频率h1=figure('numbertitle','off','name','1','pos',[50 200 420 420]);bar(t,A(t,1)),xlabel('自由度（质量块）'),ylabel('振型向量'),%输出对应的X阶主振型title('1阶主振型'),hold on,gridpause(0.1)%暂停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)%暂停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二、传递矩阵法（P45页，例2-5）function cdjzclear all %清空当前所有数据close all %关闭当前所有的绘图窗口J1=1;J2=1;J3=2; %定义各个转动惯量Jk2=1100000;k3=1200000;k4=100000; %定义各个刚度的具体数值fid=fopen('chuandi','wt'); %建立打开速度文件M1L=0; %定义变量M1L的初始值为0，M1L为第一个质量块左边的转矩大小for WN=0:.01:2000 %for 循环shita1R=1;M1R=-WN^2*J1; %定义θ1，同时利用公式2-73求出质量块右边的转矩大小M1Rshita2R=shita1R+1/k2*M1R;M2R=shita1R*(-WN^2*J2)+(1+(-WN^2*J2)/k2)*M1R;%求θ2，M2R，应用到2-78公式，问题是这个公式是θL的shita3R=shita2R+1/k2*M2R;M3R=shita2R*(-WN^2*J3)+(1+(-WN^2*J3)/k3)*M2R; %依次类推shita4R=shita3R+1/k4*M3R; %（待分析，不知神马东东）if abs(shita4R)<0.005 %abs()函数，对括号内的所有元素平方后取绝对值；判断精度是否满足，满足后执行WN %搜索到的固有频率（rad/s),并显示shita=[shita1R;shita2R;shita3R;shita4R]%搜索到阵型，并显示bar(shita),xlabel('对应的质量块'),ylabel('振型向量') %画θ矩阵图形pause(1.0) %1秒一暂停end %结束if语句fprintf(fid,'%30.15f',shita4R); %将θ4输到文档中，数据类型是30位数，15位小数的浮点数end %结束for循环fid=fopen('chuandi','rt'); % 此时fid有返回值，当是正数时代表打开文件成功，-1代表失败，rt表示读写x=fscanf(fid,'%f',[1,200001]); %要求输入1到2000001的数t=0:.01:2000;plot(t,x),grid,xlabel('频率rad/s'),ylabel('第四个质量块的转角(rad)'),title('用传递矩阵法求固有频率') 作业（P48-P49）2-8作业题function jzdd2 %第2-8作业题clear all %清空当前所有数据close all %关闭当前所有的绘图窗口fid1=fopen('A-2','wt'); %建立第一个名为“A-1”的可写文档fid2=fopen('B-2','wt'); %同上M(1,1)=1;M(2,2)=1;M(3,3)=1; %以上三段代码是为了输入质量矩阵，依题意均为1K(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; %输入刚度矩阵，依题意均为1D=inv(K)*M; %inv表示对矩阵取逆，公式2-65A=ones(3,1); %定义一个初始迭代阵型，ones()函数表示是3个位为1的单列矩阵，ones(i,j)%则是i行两列都是j都是1的数组！在这方法中一般取ones(i,1),i=质量各数for i=1:3 %进入for循环，将要执行书本上37业的迭代，循环次数为1到i，每执行一次循环代表一个质量块的主阵型pp0=0; %令初始的P（0）=0，并且输出Ｐ（０）i %输出iB=D*A; %代入2-60的公式pp=1.0/B(3); %同上，B(3),是你当时定义主阵型时，第i个元素全为1，当然也可以是第1个A=B/B(3); %同上%以上三段代码，是初始赋值while abs((pp-pp0)/pp)>1e-6 % 判断是否满足迭代的条件|P（k）^2-P(k-1)^2|<δ%满足后退出while循环pp0=pp; %P(k)^2=上诉最终的pp值B=D*A;pp=1.0/B(3);A=B/B(3);%真正地迭代代码end %结束while循环f=sqrt(pp) %公式2-61fprintf(fid1,'%20.5f',A); %fid为文件句柄，指定要写入数据的文件，format是用来控制所写数据格式的格式符，与fscanf函数相同，A是用来存放数据的矩阵。