动力学程序讲解Word格式.docx

1、%输出Aphi=%.3gn,phi）;%输出elseif z=1 %临界阻尼情况，按照1-14到1-15公式计算 a1=x0;% C1 a2=v0+wn*x0;%C2a1=%.3gn,a1）;%输出a1a2=%.3gn,a2）;%输出a2 x=（a1+a2*t）.*exp（-wn*t）;else %过阻尼，按照1-13公式计算 a1=（-v0+（-z+sqrt（z2-1）*wn*x0）/2/wn/sqrt（z2-1）; a2=（v0+（z+sqrt（z2-1）*wn*v0）/2/wn/sqrt（z2-1）; x=exp（-z*wn*t）.*（a1*exp（-wn*sqrt（z2-1）*t）+a

2、2*exp（wn*sqrt（z2-1）*t）;endplot（t,x）,grid onxlabel（时间（s）ylabel（位移（m）title（位移相对时间的关系单自由度系统的谐迫振动（P11业）function vtb2（m,c,k,x0,v0,tf,w,f0）%单自由度系统的谐迫振动lan=w/wn %的求法1-18公式 %输出Wn%输出 %输出Wda=sqrt（v0+z*wn*x0）2+（x0*wd）2）/wd2）;tf/1000:phi1=atan（x0*wd/（v0+z*wn*x0）;%按1-12求phi2=atan（2*z*lan/（1-lan2）;%求b=wn2*f0/k/sq

3、rt（wn2-w2）+（2*z*wn*w）2）;%求B的稳态响应的振幅x=a*exp（-z*wn*t）.*sin（sqrt（1-z2）*wn*t+phi1）+b*sin（w*t-phi2）;%响应方程plot（real（t）,real（x）,grid位移（cm）位移与时间的关系第二章一、矩阵迭代法（P38业，例2-3）function jzdd %矩阵迭代法clear all %清空当前所有数据close all %关闭当前所有的绘图窗口fid1=fopen（A-1,wt）; %建立第一个名为“A-1”的可写文档fid2=fopen（B-1%同上M（1,1）=2;M（2,2）=1.5;M（3,

4、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）

5、=0，并且输出（） i %输出i B=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值 %真正地迭代代码 end %结束while循环f=sqrt（pp）/2/pi %公式2-61fprintf（fid1,%20.5f %fid为文件句柄，指定要写入数据的文件，forma

6、t是用来控制所写数据格式的格式符，与fscanf函数相同，A是用来存放数据的矩阵。fprintf（fid2,f）; %储存fD=D-A*A*M/（A*M*A*pp）; %下一阶的动力矩阵【D*】，A代表A的逆矩阵end%结束for循环rt%读取A-1文档A=fscanf（fid1,%f,3,3）;%输入3x3的A矩阵，来源为fid1%读取B-1文档f=fscanf（fid2,3,1）;%输入（拾取）3x1的f矩阵,来源是fid2t=1:3;h1=figure（numbertitleoffname0pos,50 200 420 420）;bar（t,f（t,1）,xlabel（频率阶级次）,yl

7、abel（Hz）,固有频率）,hold on,grid %bar代表绘制长条图t代表起点横坐标，f（t,1）纵坐标,输出各阶频率阶次的的固有频率1bar（t,A（t,1）,xlabel（自由度（质量块）振型向量）,%输出对应的X阶主振型1阶主振型）,hold on,gridpause（0.1）%暂停0.1秒2bar（t,A（t,2）,xlabel（2阶主振型3bar（t,A（t,3）,xlabel（3阶主振型二、传递矩阵法（P45页，例2-5）function cdjzJ1=1;J2=1;J3=2; %定义各个转动惯量Jk2=1100000;k3=1200000;k4=100000; %定义各

8、个刚度的具体数值fid=fopen（chuandi %建立打开速度文件M1L=0; %定义变量M1L的初始值为0，M1L为第一个质量块左边的转矩大小for WN=0:.01:2000 %for 循环 shita1R=1;M1R=-WN2*J1; %定义1，同时利用公式2-73求出质量块右边的转矩大小M1R shita2R=shita1R+1/k2*M1R;M2R=shita1R*（-WN2*J2）+（1+（-WN2*J2）/k2）*M1R;%求2，M2R，应用到2-78公式，问题是这个公式是L的 shita3R=shita2R+1/k2*M2R;M3R=shita2R*（-WN2*J3）+（1

9、+（-WN2*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（对应的质量块） %画矩阵图形 pause（1.0） %1秒一暂停 end %结束if语句 fprintf（fid,%30.15f,shita4R）; %将4输到文档

10、中，数据类型是30位数，15位小数的浮点数end %结束 for循环 % 此时fid有返回值，当是正数时代表打开文件成功，-1代表失败，rt表示读写x=fscanf（fid,1,200001）; %要求输入1到2000001的数2000;plot（t,x）,grid,xlabel（频率rad/s第四个质量块的转角（rad）,title（用传递矩阵法求固有频率作业（P48-P49）2-8作业题function jzdd2 %第2-8作业题A-2B-2M（1,1）=1;M（2,2）=1; %以上三段代码是为了输入质量矩阵，依题意均为1K（1,1）=2;K（1,2）=-1;K（2,1）=-1;K（2

11、,2）=2; %输入刚度矩阵，依题意均为1 %定义一个初始迭代阵型，ones（）函数表示是3个位为1的单列矩阵，ones（i,j） %则是i行两列都是j都是1的数组！3 %进入for循环，将要执行书本上37业的迭代，循环次数为1到i，每执行一次循环代表一个质量块的主阵型 %同上，B（3）,是你当时定义主阵型时，第i个元素全为1，当然也可以是第1个f=sqrt（pp） %公式2-61%下一阶的动力矩阵【D*】，A代表A的逆矩阵end %结束for循环%读取A-2文档 %输入3x3的A矩阵，来源为fid1 %读取B-2文档）, %输出对应的X阶主振型）， 2-9作业题function jzdd3A

12、-3 %建立第一个名为“A-3”的可写文档B-3M（1,1）=4;M（2,2）=2;K（1,1）=4;%输入刚度矩阵 pp=1.0/B（1）; %同上，B（1）,是你当时定义主阵型时，第i个元素全为1，当然也可以是第一个 A=B/B（1）; ， 满足后退出while循环 end%结束while循环%储存f%读取A-3文档2-14作业题2-14作业题function cdjz2 %clear all, close all%清空当前所有数据，关闭当前所有的绘图窗口J1=1/2;%定义各个转动惯量J，I=1k1=100000;k2=100000;%定义各个刚度的具体数值） %建立打开速度文件M1L=

13、0 %定义变量M1L的初始值为0，M1L为第一个质量块左边的转矩大小%定义1=1，同时利用公式2-73求出质量块右边的转矩大小M1R shita2R=shita1R+1/k1*M1R;M2R=shita1R*（-WN2*J2）+（1+（-WN2*J2）/k1）*M1R; if abs（shita3R）0.005 %判断精度是否满足，满足后执行 WN %搜索到的固有频率（rad/s），并显示shita3R %搜索到振型，并显示）%画矩阵图形，条形图 pause（0.1） %1秒一暂停 endfprintf（fid,shita3R）;%储存shita3R %打开可读写chuandi文件 %拾取数据

14、200001;plot（.01*t,x）,grid,xlabel（第三个质量块的转角（rad）第三章一、欧拉法（P52业，例3-1）function vtb3（m,c,k,x0,v0,tf,w,f0,delt） %请运行vtb3（18.2,1.49,43.8,1,1,100,15,44.5,1） %用欧拉法计算单自由度系统谐迫振动响应 wn=sqrt（k/m）; %计算固有频率wn fid1=fopen（disp %建立一个位移文件disp.dat for t=0:delt:tf %delt为时间步长 xdd=（f0*sin（w*t）-k*x0-c*v0）/m; %计算加速度，代入公式3-1第

15、三条公式,自下而上 xd=v0+xdd*delt; %计算速度, x=x0+xd*delt; %计算位移x fprintf（fid1,%10.4f,x0）; %向文件中写入数据 x0=x; v0=xd; t fid2=fopen（ %打开disp.dat文件 n=tf/delt; %disp.dat文件中位移的个数 x=fscanf（fid2,1,n）; %将disp.dat文件中位移写成矩阵 t=1:n; plot（t,x）,grid xlabel（ ylabel（ title（二、欧拉法的改进（P52页）function vtb4（m,c,k,x0,v0,tf,w,f0,delt） %请运行vtb4（282,249,43.8,1,1,100,15,4.5,1）%用改进的欧拉法计算单自由度系统谐迫振动响应tf %delt